Sys.time()
rm(list = ls())
setwd([PATH])
source("0_Initiation_Programs.R",encoding = "UTF-8") ; gc()
cl<-makeCluster([CORES],type="SOCK", outfile="")
registerDoSNOW(cl)

###
#  Toggles for each set of analysis below
###
Baseline_Estimation<-T

UC_Campus_Ests<-T
STEM_Perform_Persist<-T
App_Charts_and_Discontinuities<-T
AppChanges<-T
IPUMS_Data<-T
CK_Replication<-T
Value_Added_Coefficients<-T
Miscellaneous<-T




if(Baseline_Estimation){
  load(paste0(secure_derived,"AffAct_Data.Rda"))
  
  ###
  #  Set up the data
  ###
  
  A$Asian<-A$CETHNICA%in%c("D","E","G","H","L","N","V")
  A$Hispanic<-A$CETHNICA%in%c("C","J")
  A$Black<-A$CETHNICA%in%c("B")
  A$CETHNICA_Cat<-A$CETHNICA_Cat_Pred ; A$CETHNICA_Cat[A$CETHNICA_Cat=="White"]<-"AAA" ; A$CETHNICA_Cat[A$CETHNICA_Cat%in%c("Decline","Other")]<-NA #Use predicted ethnicities for decline-to-respond applicants in some circumstances
  
  #Income quantiles by ethnicity
  A$CETHNICA_Cat_Simple<-A$CETHNICA_Cat_Pred ; A$CETHNICA_Cat_Simple[A$CETHNICA%in%c("A","B","C","J")|A$CETHNICA_Cat_Pred%in%c("Black","Hispanic")]<-"URM"
  load("Data/Derived/Income_Quants_AA.Rda")
  for(y in 2001:2017){
    for(e in c("Asian","URM","White")){
      b<-cuts[cuts$Year==y&cuts$Eth==e,]
      b$Val[1]<-0 ; b$Val[101]<-1000000000000
      b<-b[!duplicated(b$Val,fromLast=T),]
      b_fixed<-cuts[cuts$Year==2017&cuts$Eth==e,]
      b_fixed$Val[1]<-0 ; b_fixed$Val[101]<-1000000000000
      b_fixed<-b_fixed[!duplicated(b_fixed$Val,fromLast=T),]
      for(i in 6:16){
        if(!((y-i)%in%1994:2002)) next
        check<-A$CETHNICA_Cat_Simple%in%c("Other","Decline",e)&A$YEARAPAY==(y-i) #White includes other
        A[check,paste0("wage_sum_",i,"_Rank")]<-as.integer(as.character(cut(A[check,paste0("wage_sum_",i)],breaks=b$Val,labels=b$Quant[-1],include.lowest = T)))
        A[check,paste0("wage_sum_",i,"_Fixed_Rank")]<-as.integer(as.character(cut(A[check,paste0("wage_sum_",i)],breaks=b_fixed$Val,labels=b_fixed$Quant[-1],include.lowest = T)))
      }
    }
  }
  (mean(cuts$Val[cuts$Eth=="URM"&cuts$Quant==98])-mean(cuts$Val[cuts$Eth=="URM"&cuts$Quant==3]))/96 #This is the mean percentile change for URM: $1,941
  
  A<-A[A$CQUARTERAP==2&A$SCHTYPE%in%c("A","B")&!is.na(A$hsgpa),] #Restrict to main sample
  A$Prop209<-A$YEARAPAY>1997 ; A$Prop209[!inrange(A$YEARAPAY,1994,2001)]<-NA
  A$logincome_parent<-log(1+A$income_parent)
  A$VETERAN_FL<-A$VETERAN_FL=="Y"
  A$Age<-A$YEARAPAY-A$BYear
  
  A$UnivGrad<-(A$UnivGrad_NSCOnly==1)|NA_to_F(A$HASGRAD5==1)
  A$UnivGrad_SixYears<-(A$UnivGrad_SixYears_NSCOnly==1)|NA_to_F(A$HASGRAD6==1)
  A$UnivGrad_FourYears<-(A$UnivGrad_FourYears_NSCOnly==1)|NA_to_F(A$HASGRAD4==1)
  
  #Quartile indicators using the 96-97 distribution of URM AI
  A$S_Q[inrange(A$S,-10000,summary(A$S[A$YEARAPAY%in%1996:1997&A$URM])[2])]<-1
  A$S_Q[inrange(A$S,summary(A$S[A$YEARAPAY%in%1996:1997&A$URM])[2],summary(A$S[A$YEARAPAY%in%1996:1997&A$URM])[3])]<-2
  A$S_Q[inrange(A$S,summary(A$S[A$YEARAPAY%in%1996:1997&A$URM])[3],summary(A$S[A$YEARAPAY%in%1996:1997&A$URM])[5])]<-3
  A$S_Q[inrange(A$S,summary(A$S[A$YEARAPAY%in%1996:1997&A$URM])[5],10000)]<-4
  #Now do the same for Black and Hispanic separately
  A$S_Q_Black[inrange(A$S,-10000,summary(A$S[A$YEARAPAY%in%1996:1997&A$Black])[2])]<-1
  A$S_Q_Black[inrange(A$S,summary(A$S[A$YEARAPAY%in%1996:1997&A$Black])[2],summary(A$S[A$YEARAPAY%in%1996:1997&A$Black])[3])]<-2
  A$S_Q_Black[inrange(A$S,summary(A$S[A$YEARAPAY%in%1996:1997&A$Black])[3],summary(A$S[A$YEARAPAY%in%1996:1997&A$Black])[5])]<-3
  A$S_Q_Black[inrange(A$S,summary(A$S[A$YEARAPAY%in%1996:1997&A$Black])[5],10000)]<-4
  A$S_Q_Hispanic[inrange(A$S,-10000,summary(A$S[A$YEARAPAY%in%1996:1997&A$Hispanic])[2])]<-1
  A$S_Q_Hispanic[inrange(A$S,summary(A$S[A$YEARAPAY%in%1996:1997&A$Hispanic])[2],summary(A$S[A$YEARAPAY%in%1996:1997&A$Hispanic])[3])]<-2
  A$S_Q_Hispanic[inrange(A$S,summary(A$S[A$YEARAPAY%in%1996:1997&A$Hispanic])[3],summary(A$S[A$YEARAPAY%in%1996:1997&A$Hispanic])[5])]<-3
  A$S_Q_Hispanic[inrange(A$S,summary(A$S[A$YEARAPAY%in%1996:1997&A$Hispanic])[5],10000)]<-4
  
  #Wage variables
  for(v in names(A)[grepl("^wage_sum_[0-9]*$",names(A))]){
    A[,v]<-winsor(A[,v],0.01) #Winsorize 1%
    A[,paste0(v,"_Uncondit")]<-A[,v] ; A[is.na(A[,paste0(v,"_Uncondit")]),paste0(v,"_Uncondit")]<-0
    A[,paste0(v,"_Employed")]<-NA_to_F(A[,v]>0)
    A[,paste0(v,"_75000")]<-NA_to_F(A[,v]>=75000)
    A[,paste0(v,"_100000")]<-NA_to_F(A[,v]>=100000)
    A[,paste0(v,"_150000")]<-NA_to_F(A[,v]>=150000)
    A[NA_to_F(A[,v]!=0),paste0(v,"_log")]<-log(A[NA_to_F(A[,v]!=0),v])
    A[,paste0(v,"_logUncond")]<-A[,paste0(v,"_log")] ; A[is.na(A[,paste0(v,"_logUncond")]),paste0(v,"_logUncond")]<-0
  }
  A$Total_Wages_NumYears<-apply(A[,grepl("wage_sum_([6-9]|1[0-6])_Uncondit",names(A))],1,function(x) sum(x>0))
  A$Total_Wages_logUncond<-log(1+apply(A[,grepl("wage_sum_([6-9]|1[0-6])_Uncondit",names(A))],1,function(x) sum(x))/14)
  A$Total_Wages_logCondit<-log(1+apply(A[,grepl("wage_sum_([6-9]|1[0-6])_Uncondit",names(A))],1,function(x) sum(x))/A$Total_Wages_NumYears)
  A$Total_Wages_Condit<-apply(A[,grepl("wage_sum_([6-9]|1[0-6])_Uncondit",names(A))],1,function(x) sum(x))/A$Total_Wages_NumYears
  A$Total_Wages_Rank_Condit<-apply(A[,grepl("wage_sum_([6-9]|1[0-6])_Rank",names(A))],1,function(x) mean(x,na.rm=T))
  A$Total_Wages_Rank_Fixed_Condit<-apply(A[,grepl("wage_sum_([6-9]|1[0-6])_Fixed_Rank",names(A))],1,function(x) mean(x,na.rm=T))
  A$Total_Wages_Num75000<-apply(A[,grepl("wage_sum_([6-9]|1[0-6])_Uncondit",names(A))],1,function(x) sum(x>75000))
  A$Total_Wages_Num100000<-apply(A[,grepl("wage_sum_([6-9]|1[0-6])_Uncondit",names(A))],1,function(x) sum(x>100000))
  A$Total_Wages30s_NumYears<-apply(A[,grepl("wage_sum_(1[23456])_Uncondit",names(A))],1,function(x) sum(x>0))
  A$Total_Wages30s_logUncond<-log(1+apply(A[,grepl("wage_sum_(1[23456])_Uncondit",names(A))],1,function(x) sum(x))/5)
  A$Total_Wages30s_logCondit<-log(1+apply(A[,grepl("wage_sum_(1[23456])_Uncondit",names(A))],1,function(x) sum(x))/A$Total_Wages30s_NumYears)
  A$Total_Wages30s_Condit<-apply(A[,grepl("wage_sum_(1[23456])_Uncondit",names(A))],1,function(x) sum(x))/A$Total_Wages30s_NumYears
  A$Total_Wages30s_Rank_Condit<-apply(A[,grepl("wage_sum_(1[23456])_Rank",names(A))],1,function(x) mean(x,na.rm=T))
  A$Total_Wages30s_Num75000<-apply(A[,grepl("wage_sum_(1[23456])_Uncondit",names(A))],1,function(x) sum(x>75000))
  A$Total_Wages30s_Num100000<-apply(A[,grepl("wage_sum_(1[23456])_Uncondit",names(A))],1,function(x) sum(x>100000))
  
  #Define a UC-only STEM variable
  stem<-unlist(read.csv("Data/Raw/Definition of STEM.csv",stringsAsFactors = F)) ; stem<-gsub("^.*[^0-9]([0-9]{6})[^0-9].*$","\\1",gsub("[.]","",stem[grepl("[.]",stem)])) ; stem[nchar(stem)==5]<-paste0("0",stem[nchar(stem)==5]) ;  for(z in stem[grepl("X",stem)]) stem<-c(stem,paste0(gsub("X","",z),addzeros(1:9999,4)))
    A$STEM_NoNSC<-NA_to_F(A$CIP6RECENT%in%stem) ; A$STEM_NoNSC[is.na(A$HASGRAD6)]<-NA
    A$STEM_NoNSC_Condit<-A$STEM_NoNSC ; A$STEM_NoNSC_Condit[NA_to_T(A$HASGRAD6==0)]<-NA
    A$STEM_UCNSC<-NA_to_F(A$STEM_NoNSC==1)|NA_to_F(A$STEM_Uncondit_NSCOnly==1)
    A$STEM_UCNSC_Condit<-A$STEM_UCNSC ; A$STEM_UCNSC_Condit[NA_to_T(A$UnivGrad_SixYears==0)]<-NA
  
  #Define Major variables
  for(v in names(table(A$NSC_Grad_MajorSimp[A$YEARAPAY%in%1996:1999])[table(A$NSC_Grad_MajorSimp[A$YEARAPAY%in%1996:1999])>1000])){
    A[,paste0("NSC_Major_",gsub("[ ']","",v))]<-NA_to_F(A$NSC_MajorSimp1==v)|NA_to_F(A$NSC_MajorSimp2==v)|NA_to_F(A$NSC_MajorSimp3==v)
    A[,paste0("NSC_Major_",gsub("[ ']","",v),"_Condit")]<-A[,paste0("NSC_Major_",gsub("[ ']","",v))] ; A[A$UnivGrad_SixYears_NSCOnly==0,paste0("NSC_Major_",gsub("[ ']","",v),"_Condit")]<-NA
  }
  
  #Derive 1-to-1 comparable VA and wage variables
  for(i in c("Mountjoy","DKCB","Chetty")){
    A[,paste0("UnivGrad_SixYears_",i)]<-A$UnivGrad_SixYears ; A[is.na(A[,paste0("FE_GradVA_",i,"_AA")]),paste0("UnivGrad_SixYears_",i)]<-NA
    A[,paste0("Total_Wages30s_Condit_",i)]<-A$Total_Wages30s_Condit ; A[is.na(A[,paste0("FE_WageVA_",i,"_AA")]),paste0("Total_Wages30s_Condit_",i)]<-NA
    A[,paste0("FE_WageVA_",i,"_AA_Condit")]<-A[,paste0("FE_WageVA_",i,"_AA")] ; A[is.na(A$Total_Wages30s_Condit),paste0("FE_WageVA_",i,"_AA_Condit")]<-NA
    for(j in c("_Eth")){
      A[,paste0("UnivGrad_SixYears_",i,j)]<-A$UnivGrad_SixYears ; A[is.na(A[,paste0("FE_GradVA_",i,"_AA",j)]),paste0("UnivGrad_SixYears_",i,j)]<-NA
      A[,paste0("Total_Wages30s_Condit_",i,j)]<-A$Total_Wages30s_Condit ; A[is.na(A[,paste0("FE_WageVA_",i,"_AA",j)]),paste0("Total_Wages30s_Condit_",i,j)]<-NA
      A[,paste0("FE_WageVA_",i,"_AA",j,"_Condit")]<-A[,paste0("FE_WageVA_",i,"_AA",j)] ; A[is.na(A$Total_Wages30s_Condit),paste0("FE_WageVA_",i,"_AA",j,"_Condit")]<-NA
    }
  }
  
  #Add in number of annual apps by HS-by-Eth-by-AIS in 94/95
  A$AIS_rounded<-round(A$AIS/500)*500 ; A$hsgpa_rounded<-winsor(round(A$hsgpa/.25)*.25,.02) ; A$Count<-1 ; A$YEARAPAY_rounded<-A$YEARAPAY-(A$YEARAPAY==1995) #-(A$YEARAPAY%%2==1)
  A<-A[,!grepl("NumApps_HS_(AIS|hsgpa)",names(A))]
  hi<-aggregate(Count~CPREVSCH+EthGen+AIS_rounded+YEARAPAY,A[A$NotIneligible&A$CQUARTERAP==2&A$adm_type=="Freshman"&A$SEX%in%c("F","M")&!grepl("Decline|Other",A$EthGen),],sum) ; names(hi)[5]<-"NumApps_HS_AIS" ; A<-merge(A,hi,all.x=T)
  hi<-aggregate(Count~CPREVSCH+EthGen+AIS_rounded,A[A$YEARAPAY%in%1994:1995&A$NotIneligible&A$CQUARTERAP==2&A$adm_type=="Freshman"&A$SEX%in%c("F","M")&!grepl("Decline|Other",A$EthGen),],sum) ; hi$Count<-hi$Count/2 ; names(hi)[4]<-"NumApps_HS_AIS9495" ; A<-merge(A,hi,all.x=T)
  hi<-aggregate(Count~CPREVSCH+EthGen+hsgpa_rounded+YEARAPAY,A[A$NotIneligible&A$CQUARTERAP==2&A$adm_type=="Freshman"&A$SEX%in%c("F","M")&!grepl("Decline|Other",A$EthGen),],sum) ; names(hi)[5]<-"NumApps_HS_hsgpa" ; A<-merge(A,hi,all.x=T)
  hi<-aggregate(Count~CPREVSCH+EthGen+hsgpa_rounded,A[A$YEARAPAY%in%1994:1995&A$NotIneligible&A$CQUARTERAP==2&A$adm_type=="Freshman"&A$SEX%in%c("F","M")&!grepl("Decline|Other",A$EthGen),],sum) ; hi$Count<-hi$Count/2 ; names(hi)[4]<-"NumApps_HS_hsgpa9495" ; A<-merge(A,hi,all.x=T)
  
  #Define control functions
  A$Prop_App_HS_AIS<-A$NumApps_HS_AIS/A$NumGrads_HS
  A$Prop_App_HS_AIS[is.infinite(A$Prop_App_HS_AIS)]<-NA ; A$Prop_App_HS_AIS[NA_to_F(A$Prop_App_HS_AIS>1)]<-1
  A$Prop_App_HS_LOO<-(A$NumApps_HS-1)/(A$NumGrads_HS-1)
  A$Prop_App_HS_LOO[A$NumGrads_HS<2|A$NumApps_HS<1]<-NA ; A$Prop_App_HS_LOO[NA_to_F(A$Prop_App_HS_LOO>1)]<-1
  temp<-qnorm(A$Prop_App_HS_AIS) ; temp[is.infinite(temp)]<-NA
  A$Heckit_AIS<-dnorm(temp)/pnorm(temp)
  temp<-qnorm(A$Prop_App_HS_LOO) ; temp[is.infinite(temp)]<-NA
  A$Heckit_LOO<-dnorm(temp)/pnorm(temp)
  A$Heckit_AIS_zeros<-A$Heckit_AIS ; A$Heckit_AIS_zeros[is.na(A$Heckit_AIS)]<-0 ; A$Heckit_AIS_missing<-is.na(A$Heckit_AIS)
  A$Heckit_LOO_zeros<-A$Heckit_LOO ; A$Heckit_LOO_zeros[is.na(A$Heckit_LOO)]<-0 ; A$Heckit_LOO_missing<-is.na(A$Heckit_LOO)
  
  ###
  #  Estimation begins here
  ###
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces the statistics presented in Tables II, III, IV, E-1, I-4, A-7, A-8, A-10, A-11, A-12, and A-14.
  ##########################################################################################################
  ##########################################################################################################
  #List of conditions for diff-in-diff regressions
  conds<-list(All=rep(T,nrow(A)), #Baseline estimates: diff-in-diff comparing 96-97 to 98-99 URM and non-URM applicants
              Q1=NA_to_F(A$S_Q==1), #Restrict sample to first AI quartile (defined among URM applicants)
              Q2=NA_to_F(A$S_Q==2), #Restrict sample to second AI quartile (defined among URM applicants)
              Q3=NA_to_F(A$S_Q==3), #Restrict sample to third AI quartile (defined among URM applicants)
              Q4=NA_to_F(A$S_Q==4), #Restrict sample to fourth AI quartile (defined among URM applicants)
              Asian=rep(T,nrow(A)), #Drop URM applicants and compare Asian to white
              Hispanic=rep(T,nrow(A)), #Drop Non-Hispanic URM applicants and compare to non-URM
              Black=rep(T,nrow(A)), #Drop Non-Black URM applicants and compare to URM
              Use95=rep(T,nrow(A))) #Replace 96-97 applicants with 94-95 applicants
             
  #List of conditions for URM-vs-Asian-vs-white diff-in-diff regressions and for Black-vs-Hispanic-vs-nonURM regressions
  conds_multieth_URM<-list(All=rep(T,nrow(A))) #Replace 96-97 applicants with 94-95 applicants
  
  #Changes in applicant characteristics (and see which variables may be driving the effect...)
  varsets<-list(c("Prop_App_HS_AIS","Prop_App_HS_LOO","Heckit_AIS","Heckit_LOO"),
                c("HighUCEnrFirstIncCC","MidUCEnrFirstIncCC","LowUCEnrFirstIncCC","CSUEnrFirstIncCC","CCEnrFirstIncCC","IvyEnrFirstIncCC","CAPrivEnrFirstIncCC","NonCAEnrFirstIncCC","NonEnrFirstIncCC"),
                c("UnivGrad_FourYears_NSCOnly","UnivGrad_SixYears_NSCOnly","STEM_Uncondit_NSCOnly","STEM_NSCOnly","HASGRAD4","HASGRAD6","STEM_NoNSC","STEM_NoNSC_Condit","UnivGrad_FourYears","UnivGrad_SixYears","STEM_UCNSC","STEM_UCNSC_Condit","MA","JD","MA_STEM","MD","MBA",names(A)[grepl("Total_Wages",names(A))]),
                names(A)[grepl("wage_sum_[0-9]*(_?[0-9]*|_Uncondit|_Employed|_log|_logUncond|_Rank)?$",names(A))],
                c("NSC_EnrFirst_IPEDS_AdmRate","NSC_EnrFirst_IPEDS_SAT1600","NSC_EnrFirst_IPEDS_GradRateSix","FE_GradVA_Mountjoy_AA","FE_WageVA_Mountjoy_AA","FE_GradVA_DKCB_AA","FE_WageVA_DKCB_AA","FE_GradVA_Chetty_AA","FE_WageVA_Chetty_AA","FE_GradVA_Chetty_AA_Eth","FE_WageVA_Chetty_AA_Eth","FE_GradVA_Mountjoy_AA","UnivGrad_SixYears_Mountjoy","FE_WageVA_Mountjoy_AA_Condit","Total_Wages30s_Condit_Mountjoy","FE_GradVA_DKCB_AA","UnivGrad_SixYears_DKCB","FE_WageVA_DKCB_AA_Condit","Total_Wages30s_Condit_DKCB","FE_GradVA_Chetty_AA","UnivGrad_SixYears_Chetty","FE_WageVA_Chetty_AA_Condit","Total_Wages30s_Condit_Chetty","FE_GradVA_Chetty_AA_Eth","UnivGrad_SixYears_Chetty_Eth","FE_WageVA_Chetty_AA_Eth_Condit","Total_Wages30s_Condit_Chetty_Eth"),
                names(A)[grepl("^NSC_Major_",names(A))])
  varnames<-c("Placebo","Enr","Outcomes","Wages","IPEDS","Majors")
  #Now choose which models to run
  version_run<-c(1,2,3,4) #There are four versions:
    #1: Baseline
    #2: Separate estimates for URM and Asian applicants, relative to white
    #3: Single-difference regressions for URM and non-URM
    #4: Separate estimates for Black and Hispanic, relative to non-URM
  var_run<-c(1,2,3,4,5,6) #There are six sets of variables, listed above; choose which sets to estimate
  
  for(version in version_run){
    if(!version%in%version_run) next
    multieth<-(version==2) ; singlediff<-(version==3) ; multieth_URM<-(version==4)
    m<-"" ; if(multieth) m<-"_Multieth" ; if(singlediff) m<-paste0(m,"_SingleDiff") ; if(multieth_URM) m<-paste0(m,"_Multieth_URM")
    if(length(var_run)<length(varsets)) load(paste0("Figures/Tables",m,".Rda")) #If not running everything, then load the ones that aren't being run
    for(v in 1:length(varsets)){
      if(!v%in%var_run) next
      CONDS<-conds ; if(multieth_URM) CONDS<-conds_multieth_URM
      R<-data.frame(matrix("",nrow=11*length(CONDS)+1,ncol=(length(varsets[[v]])+1)*(1+(2*(multieth|multieth_URM)))),stringsAsFactors = F)
      r<-1 ; for(z in 1:length(CONDS)){
        c<-2
        R[r,1]<-names(CONDS)[z] ; print(names(CONDS)[z])
        for(i in varsets[[v]]){
          numdec=2 ; if(v%in%c(2,5)) numdec<-1 ; if(v==4) numdec<-5
          if(grepl("Total.*Condit|^FE_Wage",i)&!grepl("Rank",i)) numdec<-0 ; if(grepl("Total.*log",i)) numdec<-3
          if((i=="Applied_Fresh2_satact")){ ; c<-c+1+(2*(multieth|multieth_URM)) ; next ; }
          if(grepl("95|96",names(CONDS)[z])){
            if(grepl("95",names(CONDS)[z])) years<-c(1994:1995,1998:1999) else years<-c(1994:1995,1996:1997)
            if(grepl("NSCOnly|UnivGrad|EnrFirst|UCNSC|UnivGrad(_SixYears)$|^NSC_Major",i)) years<-years[!years==1994]
          } else years<-1996:1999
          R<-genmat(R,r,c,i,A[A$YEARAPAY%in%(years)&CONDS[[z]],],asian=grepl("Asian",names(CONDS)[z]),hispanic=grepl("Hispanic",names(CONDS)[z]),black=grepl("Black",names(CONDS)[z]),numdec=numdec,multieth=multieth,singlediff=singlediff,multieth_URM=multieth_URM,Use9697=grepl("96",names(CONDS)[z]))
          names(R)[c]<-i ; c<-c+1+(2*(multieth|multieth_URM))
          print(i)
        }
        r<-r+11
      }
      assign(paste0("R",v-1),R) ; gc()
      write.csv(R,file=paste0("Figures/Prop209_",varnames[v],m,".csv"))
    }
    save(R0,R1,R2,R3,R4,R5,file=paste0("Figures/Tables",m,".Rda"))
  }
  
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces Figures V, E-1, A-9, and A-12.
  ##########################################################################################################
  ##########################################################################################################
  for(singlediff in c(T,F)){
  for(multieth in c("","_Multieth","_Multieth_URM")){
    if(singlediff&multieth!="") next
    CONDS<-conds ; if(multieth=="_Multieth_URM") CONDS<-conds_multieth_URM
    n<-c("Employment Rate","Percentage Points","Percentage Points","Percentage Points","Percentage Points","Dollars","Dollars","Log Dollars","Log Dollars","Percentiles")
    m<-multieth ; if(singlediff) m<-paste0(m,"_SingleDiff")
    load(paste0("Figures/Tables",m,".Rda"))
    for(ss in 1:length(CONDS)){
      i<-1
      r<-which(R3[,1]==names(CONDS)[ss])
      for(t in c("Employed","75000","100000","125000","150000","Uncondit","wage_sum","log","logUncond","Rank")){
        for(v in c("OLS")){ #,"IV"
          if(ss==2&v=="IV") next
          if(v=="OLS") rplus<-3 else rplus<-10
          d1<-t(apply(R3[(rplus:(rplus+1))+r,grepl(paste0(t,"(_[0-9]*)?$"),names(R3))],2,function(x) as.numeric(gsub("[^0-9.-]","",as.matrix(x)))))
          d1<-cbind(as.integer(gsub("^.*_([0-9]{1,2})(_.*)?$","\\1",row.names(d1))),d1)
          d<-d1[!grepl("^(Asian|URM)",row.names(d1))&!grepl("1[789]",row.names(d1)),]
          if(nrow(d)==0) next
          ylim=c(round(min(d[,2]-2*d[,3])*1.1/1000)*1000,round(max(d[,2]+2*d[,3])*1.1/1000)*1000) ; yaxp<-NA
          if(t%in%c("75000","100000","125000","150000","Employed")) ylim=c(round((min(d[,2]-2*d[,3])-1)/1)*1,round((max(d[,2]+2*d[,3])+1)/1)*1)
          if(t%in%c("log","logUncond")) ylim=c(round((min(d[,2]-2*d[,3])-.05)/.05)*.05,round((max(d[,2]+2*d[,3])+.05)/.05)*.05)
          if(t=="log"&ss%in%c(1)){ ylim<-c(-.1,.03) ; yaxp<-c(-.1,-.05,0) ; }
          if(t=="log"&ss%in%c(8)) ylim<-c(-.15,.05) #; yaxp<-c()
          if(t=="log"&ss%in%c(2,3,4,5)) ylim<-c(-.25,.15)
          if(t=="log"&ss%in%c(7)) ylim<-c(-.1,.1)
          if(t=="wage_sum"&ss%in%c(1,8)){ ylim<-c(-5000,1000) ; yaxp<-as.integer(c(-5000,-2500,0)) ; }
          if(t=="wage_sum"&ss%in%c(2,3,4,5)){ ylim<-c(-10000,5000) ; yaxp<-c(-10000,0,5000) ; }
          if(t%in%c("75000","100000","125000","150000")&ss==1) ylim<-c(-4,1)
          if(t%in%c("75000","100000","125000","150000")&ss==7) ylim<-c(-4,1.5)
          if(t%in%c("75000","100000","125000","150000")&ss==7) ylim<-c(-4,1.5)
          if(t%in%c("Rank")) ylim<-c(-5,2)
          tiff(paste0("Figures/Prop209_",t,"_",v,names(CONDS)[ss],m,".tiff"),5,3.5,units="in",res=300)
            if(is.na(yaxp)) plot(d[,1], d[,2],ylim=ylim,cex.lab=1.5-0.25*(t=="Rank"),cex.axis=1.6-0.1*(t=="Rank"),pch=19, xlab="Years After Application", ylab=paste0(n[i]),mgp=c(2.5+0.5*(t=="Rank"),1,0))
            if(!is.na(yaxp)){
              plot(d[,1], d[,2],ylim=ylim,cex.lab=1.5,cex.axis=1.6,pch=19, xlab="Years After Application", ylab=paste0(n[i]), axes=F,mgp=c(2.5,1,0))
              axis(1, labels = seq(6,16,2), at = seq(6,16,2),cex.axis=1.6)
              if(sum(inrange(abs(yaxp),0.001,0.999))==0) axis(2, labels = ex(yaxp,0,justnum = T,comma=T), at = yaxp,cex.axis=1.5) else axis(2, labels = yaxp, at = yaxp,cex.axis=1.2)
              box()
            }
            lines(c(-20,20),c(0,0),col="gray")
            arrows(d[,1], d[,2]-d[,3]*1.96, d[,1], d[,2]+d[,3]*1.96, length=0.05, angle=90, code=3)
          dev.off()
          
          #Replicate, including '94-95
          if(ss==1){
            ylim2=ylim
            if(t%in%c("75000","100000","125000","150000")&ss==1) ylim2<-c(-7,1)
            if(t=="wage_sum"&ss%in%c(1,8)){ ylim<-c(-7000,1000) ; yaxp<-as.integer(c(-5000,-2500,0)) ; }
            d1<-t(apply(R3[(rplus:(rplus+1))+78,grepl(paste0(t,"(_[0-9]*)?$"),names(R3))],2,function(x) as.numeric(gsub("[^0-9.-]","",as.matrix(x)))))
            d1<-cbind(as.integer(gsub("^.*_([0-9]{1,2})(_.*)?$","\\1",row.names(d1))),d1)
            d2<-d1[!grepl("^(Asian|URM)",row.names(d1))&!grepl("1[789]",row.names(d1)),]
            shift<-.1
            tiff(paste0("Figures/Prop209_Inc9495_",t,"_",v,names(CONDS)[ss],m,".tiff"),5,3.5,units="in",res=300)
              if(is.na(yaxp)){
                plot(d[,1]-shift, d[,2],ylim=ylim2,cex.lab=1.5,cex.axis=1.6,pch=19, xlab="Years After Application", ylab=paste0(n[i]),mgp=c(2.5,1,0))
              } 
              if(!is.na(yaxp)){
                plot(d[,1]-shift, d[,2],ylim=ylim,cex.lab=1.5,cex.axis=1.6,pch=19, xlab="Years After Application", ylab=paste0(n[i]), axes=F,mgp=c(2.5,1,0))
                axis(1, labels = seq(6,16,2), at = seq(6,16,2),cex.axis=1.6)
                if(sum(inrange(abs(yaxp),0.001,0.999))==0) axis(2, labels = ex(yaxp,0,justnum = T,comma=T), at = yaxp,cex.axis=1.5) else axis(2, labels = yaxp, at = yaxp,cex.axis=1.2)
                box()
              }
              lines(c(-20,20),c(0,0),col="gray")
              arrows(d[,1]-shift, d[,2]-d[,3]*1.96, d[,1]-shift, d[,2]+d[,3]*1.96, length=0.05, angle=90, code=3)
              points(d2[,1]+shift, d2[,2],pch=19,col="gray50")
              arrows(d2[,1]+shift, d2[,2]-d2[,3]*1.96, d2[,1]+shift, d2[,2]+d2[,3]*1.96, length=0.05, angle=90, code=3,col="gray50")
            dev.off()
          }
          
          if(multieth!=""){
            if(names(CONDS)[ss]=="Asian") next
            d1<-t(apply(R3[(rplus:(rplus+1))+r,grepl(paste0(t,"(_[0-9]*)?$"),names(R3))],2,function(x) as.numeric(gsub("[^0-9.-]","",as.matrix(x))))) #Replicated from above
            d1<-cbind(as.integer(gsub("^.*_([0-9]{1,2})(_.*)?$","\\1",row.names(d1))),d1)
            d<-d1[grepl("^(Asian|URM|Hispanic|Black)",row.names(d1))&!grepl("1[789]",row.names(d1)),]
            urm<-grepl("^URM|Black",row.names(d)) ; asian<-grepl("^Asian|Hispanic",row.names(d))
            ylim=c(round(min(d[,2]-2*d[,3])*1.1/1000)*1000,round(max(d[,2]+2*d[,3])*1.1/1000)*1000)
            if(t%in%c("75000","100000","125000","150000","Employed")) ylim=c(round((min(d[,2]-2*d[,3])-1)/1)*1,round((max(d[,2]+2*d[,3])+1)/1)*1)
            if(t%in%c("log","logUncond")) ylim=c(round((min(d[,2]-2*d[,3])-.1)/.1)*.1,round((max(d[,2]+2*d[,3])+.1)/.1)*.1)
            if(t=="log"&ss%in%c(1)) ylim<-c(-.1,.08) #; yaxp<-c()
            if(t%in%c("75000","100000","125000","150000")&ss==1) ylim<-c(-5,2)
            tiff(paste0("Figures/Prop209_",t,"_",v,names(CONDS)[ss],m,".tiff"),5,3.5,units="in",res=300)
              plot(d[urm,1]-.1, d[urm,2],
                   ylim=ylim,cex.lab=1.5,cex.axis=1.75,
                   pch=19-2*(sum(grepl("Black",row.names(d)))>0), xlab="Years After Application", ylab=paste0(n[i]),mgp=c(2.4,1,0)
              )
              points(d[asian,1]+.1, d[asian,2],
                     ylim=ylim,pch=19-4*(sum(grepl("Black",row.names(d)))>0),col="gray50")
              lines(c(-20,20),c(0,0),col="gray70")
              arrows(d[urm,1]-.1, d[urm,2]-d[urm,3]*1.96, d[urm,1]-.1, d[urm,2]+d[urm,3]*1.96, length=0.05, angle=90, code=3)
              arrows(d[asian,1]+.1, d[asian,2]-d[asian,3]*1.96, d[asian,1]+.1, d[asian,2]+d[asian,3]*1.96, length=0.05, angle=90, code=3,col="gray50")
            dev.off()
          }
        }
        i<-i+1
      }
    }
  }
  }
  #Produce legends for figures
  tiff(paste0("Figures/Prop209_9495Legend.tiff"),width=8,height=3,units="in",res=500)
    plot(100000,100000,bty="n",box=F,yaxt="n",xaxt="n",xlim=c(0,50),ylim=c(0,10),xlab="",ylab="")
    legend(1,8,legend=c("'96-97 Baseline","'94-95 Baseline"),pch=c(18,18),lty=c(1,1),lwd=c(1,1),col=c("black","gray50"),cex=.6,ncol = 2)
  dev.off()
  tiff(paste0("Figures/Prop209_AsianLegend.tiff"),width=8,height=3,units="in",res=500)
    plot(100000,100000,bty="n",box=F,yaxt="n",xaxt="n",xlim=c(0,50),ylim=c(0,10),xlab="",ylab="")
    legend(1,8,legend=c("URM","Asian"),pch=c(18,18),lty=c(1,1),lwd=c(1,1),col=c("black","gray50"),cex=.6,ncol = 2)
  dev.off()
  tiff(paste0("Figures/Prop209_MultiEth_URMLegend.tiff"),width=8,height=3,units="in",res=500)
    plot(100000,100000,bty="n",box=F,yaxt="n",xaxt="n",xlim=c(0,50),ylim=c(0,10),xlab="",ylab="")
    legend(1,8,legend=c("Black","Hispanic"),pch=c(18,18),lty=c(1,1),lwd=c(1,1),col=c("black","gray50"),cex=.6,ncol = 2)
  dev.off()
  tiff(paste0("Figures/Prop209_NonURMLegend.tiff"),width=8,height=3,units="in",res=500)
    plot(100000,100000,bty="n",box=F,yaxt="n",xaxt="n",xlim=c(0,50),ylim=c(0,10),xlab="",ylab="")
    legend(1,8,legend=c("URM","Non-URM"),pch=c(18,18),lty=c(1,1),lwd=c(1,1),col=c("black","gray50"),cex=.6,ncol = 2)
  dev.off()
  tiff(paste0("Figures/Prop209_MultiURMLegend.tiff"),width=8,height=3,units="in",res=500)
    plot(100000,100000,bty="n",box=F,yaxt="n",xaxt="n",xlim=c(0,50),ylim=c(0,10),xlab="",ylab="")
    legend(1,8,legend=c("Black","Hispanic"),pch=c(17,15),lty=c(1,1),lwd=c(1,1),col=c("black","gray50"),cex=.6,ncol = 2)
  dev.off()
  
  
  
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces Figure A-15.
  ##########################################################################################################
  ##########################################################################################################
  A$SAT_Sum<-A$SATIV+A$SATIM+A$SATII_W+A$SATII_M+A$SATII_Other
  for(v in c("UnivGrad_SixYears","STEM_UCNSC","MA","Total_Wages_logCondit","Total_Wages30s_Num100000","Total_Wages_Rank_Condit")){
    R_grad<-matrix(NA,nrow=1000,ncol=60) ; c<-1
    controls<-c("logincome_parent_zeros","HSCRSE_TYP_ACADHONRS_12","PAR_OCCUPN_KEY","Female","educ_lvl_max","Heckit_AIS_zeros","Heckit_LOO_zeros","NotIneligible","hsgpa_rank")
    data<-A[A$YEARAPAY%in%(1996:1999),]
    data$DepVar<-data[,v] ; if(grepl("UnivGrad|MA|STEM",v)) data$DepVar<-data$DepVar*100
    if(v%in%c("UnivGrad_SixYears")){ ; data$DepVar[NA_to_T(data$S_Q!=1)]<-NA ; }
    if(v%in%c("Heckit_LOO","Prop_App_HS_LOO")) controls<-controls[!grepl("Heckit",controls)]
    for(i in 1:length(controls)){
      hi<-combn(controls,i)
      r<-1
      for(j in 1:ncol(hi)){
        for(zz in 1:nrow(hi)) data[,paste0("Int",zz)]<-data[,hi[zz,j]]
        if("logincome_parent_zeros"%in%hi[,j]){ ; zz<-zz+1 ; data[,paste0("Int",zz)]<-data$income_Missing ; }
        if("hsgpa_rank"%in%hi[,j]){ ; zz<-zz+1 ; data[,paste0("Int",zz)]<-data$APPL_SLFRPTED_FRSHMN_GPA_SCO ; }
        if("SATIV"%in%hi[,j]){ ; zz<-zz+1 ; data[,paste0("Int",zz)]<-data$SATIM ; }
        satii<-"" ; if("SATIIW"%in%hi[,j]){ ; zz<-zz+3 ; data[,paste0("Int",zz-2)]<-data$SATII_M ; data[,paste0("Int",zz-1)]<-data$SATII_M2_Indic ; data[,paste0("Int",zz)]<-data$SATII_Other ; satii<-"+SATII_1" ; }
        if("Heckit_LOO_zeros"%in%hi[,j]){ ; zz<-zz+1 ; data[,paste0("Int",zz)]<-data$Heckit_LOO_missing ; }
        if("Heckit_AIS_zeros"%in%hi[,j]){ ; zz<-zz+1 ; data[,paste0("Int",zz)]<-data$Heckit_AIS_missing ; }
        occ<-"" ; if("PAR_OCCUPN_KEY"%in%hi[,j]){ ; data[,paste0("Int",which(hi[,j]=="PAR_OCCUPN_KEY"))]<-1 ; occ<-"+PAR_OCCUPN_KEY" ; }
        ed<-"" ; if("educ_lvl_max"%in%hi[,j]){ ; data[,paste0("Int",which(hi[,j]=="educ_lvl_max"))]<-1 ; ed<-"+educ_lvl_max" ; }
        
        t<-summary(felm(formula(paste0("DepVar~URM*Prop209+SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m+",paste0("Int",1:zz,collapse="+"),"| factor(CPREVSCH)+factor(SATII_1_nom)",occ,ed,satii,"")),data=data,exactDOF = T),robust=T) #,hs
        R_grad[r,c:(c+2)]<-t$coefficients[which(row.names(t$coefficients)=="URMTRUE:Prop209TRUE"),1:3]
        R_grad[r,c+3]<-t$r2
        r<-r+1
      }
      print(paste(Sys.time(),i))
      c<-c+4
    }
    r_grad<-data.frame(Column=1:(length(controls)+3),beta_min=rep(NA,(length(controls)+3)),sterr_min=rep(NA,(length(controls)+3)),t_min=rep(NA,(length(controls)+3)),beta_max=rep(NA,(length(controls)+3)),sterr_max=rep(NA,(length(controls)+3)),t_max=rep(NA,(length(controls)+3)),beta_median=rep(NA,(length(controls)+3)),r_min=rep(NA,(length(controls)+3)),r_max=rep(NA,(length(controls)+3)),r_median=rep(NA,(length(controls)+3)))
    t<-summary(felm(DepVar~URM*Prop209,data=data),robust=T)
      r_grad[1,2:7]<-rep(t$coefficients[which(row.names(t$coefficients)=="URMTRUE:Prop209TRUE"),1:3],2)
      r_grad[1,9:11]<-rep(t$r2,3)
    t<-summary(felm(DepVar~URM*Prop209 | factor(CPREVSCH),data=data),robust=T)
      r_grad[2,2:7]<-rep(t$coefficients[which(row.names(t$coefficients)=="URMTRUE:Prop209TRUE"),1:3],2)
      r_grad[2,9:11]<-rep(t$r2,3)
    t<-summary(felm(DepVar~URM*Prop209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom),data=data),robust=T)
      r_grad[3,2:7]<-rep(t$coefficients[which(row.names(t$coefficients)=="URMTRUE:Prop209TRUE"),1:3],2)
      r_grad[3,9:11]<-rep(t$r2,3)
    r_grad[4:(length(controls)+3),2]<-apply(R_grad[,seq(1,length(controls)*4,4)],2,function(x) min(x,na.rm=T))
    r_grad[4:(length(controls)+3),4]<-apply(R_grad[,seq(3,length(controls)*4,4)],2,function(x) min(x,na.rm=T))
    r_grad[4:(length(controls)+3),5]<-apply(R_grad[,seq(1,length(controls)*4,4)],2,function(x) max(x,na.rm=T))
    r_grad[4:(length(controls)+3),7]<-apply(R_grad[,seq(3,length(controls)*4,4)],2,function(x) max(x,na.rm=T))
    r_grad[4:(length(controls)+3),8]<-apply(R_grad[,seq(1,length(controls)*4,4)],2,function(x) median(x,na.rm=T)) ; r_grad[3,8]<-r_grad[3,2]
    r_grad[4:(length(controls)+3),9]<-apply(R_grad[,seq(4,length(controls)*4,4)],2,function(x) min(x,na.rm=T))
    r_grad[4:(length(controls)+3),10]<-apply(R_grad[,seq(4,length(controls)*4,4)],2,function(x) max(x,na.rm=T))
    r_grad[4:(length(controls)+3),11]<-apply(R_grad[,seq(4,length(controls)*4,4)],2,function(x) median(x,na.rm=T))
    for(i in 1:length(controls)){
      r_grad[i+3,3]<-R_grad[which(R_grad[,i*4-3]-r_grad[i+3,2]<0.000001)[1],i*4-2]
      r_grad[i+3,6]<-R_grad[which(R_grad[,i*4-3]-r_grad[i+3,5]<0.000001)[1],i*4-2]
    }
    r_grad<-r_grad[-1,] ; r_grad$Column<-r_grad$Column-2
    ylabel<-"# of Years" ; yaxl<-NA
    if(grepl("Total_Wages_Num75000",v)){ ; ylim<-c(-.36,.1) ; y_ax<-round(seq(-.3,.1,.1),1)
    } else if(grepl("Total_Wages_Num100000",v)){ ; ylim<-c(-.36,.1) ; y_ax<-round(seq(-.3,.1,.1),1)
    } else if(grepl("Total_Wages_logCondit",v)){ ; ylim<-c(-.075,.01) ; y_ax<-round(seq(-.075,0,.025),3) ;  yaxl<-c("","-0.05","","0") ; ylabel<-"Log Dollars"
    } else if(grepl("Total_Wages30s_Num75000",v)){ ; ylim<-c(-.15,.05) ; y_ax<-round(seq(-.15,.05,.05),2)
    } else if(grepl("Total_Wages30s_Num100000",v)){ ; ylim<-c(-.11,.01) ; y_ax<-round(seq(-.1,0,.025),3) ; yaxl<-c("-0.1","","-0.05","","0")
    } else if(grepl("Total_Wages30s_logCondit",v)){ ; ylim<-c(-.075,.025) ; y_ax<-round(seq(-.075,.025,.025),3) ; ylabel<-"Log Dollars"
    } else if(grepl("Total_Wages_Rank_Condit",v)){ ; ylim<-c(-2,.1) ; y_ax<-round(seq(-2,0,1),0) ; ylabel<-"Percentiles"
    } else if(grepl("Total_Wages_Condit",v)){ ; ylim<-c(-.075,.01) ; y_ax<-round(seq(-.075,0,.025),3) ;  yaxl<-c("","-0.05","","0") ; ylabel<-"Log Dollars"
    } else if(grepl("Total_Wages30s_Condit",v)){ ; ylim<-c(-4000,500) ; y_ax<-round(seq(-4000,2000,2000),0) ; ylabel<-"Dollars"
    } else if(grepl("STEM",v)){ ; ylim<-c(-2,1) ; y_ax<-round(seq(-2,1,1),1) ; ylabel<-"Percent"
    } else if(grepl("MA",v)){ ; ylim<-c(-3,1) ; y_ax<-seq(-3,1,1) ; ylabel<-"Percent"
    } else if(grepl("Heckit",v)){ ; ylim<-c(-0.1,0.1) ; y_ax<-seq(-0.1,0.1,.05) ; ylabel<-"Units"
    } else if(grepl("Prop_App_HS_LOO",v)){ ; ylim<-c(-0.05,0.02) ; y_ax<-seq(-0.05,0.02,.01) ; ylabel<-"Percent"
    } else{ ; ylim<-c(-9,0.5) ; y_ax<-seq(-8,2,2) ; ylabel<-"Percent" ; }
    tiff(paste0("Figures/Prop209_",v,"_Robustness.tiff"),5,3.5,"in",res=300) #,width = 1000,height=1000
      plot(r_grad$Column,r_grad$beta_min,
           ylim=ylim,
           pch=19, xlab="Specification: Number of Covariates", ylab=ylabel,axes=F,cex.lab=1.25,mgp=c(2.5,1,0))
      lines(c(-20,20),c(0,0),col="gray")
      lines(r_grad$Column,r_grad$beta_min,col="gray40")
      points(r_grad$Column,r_grad$beta_max,pch=19)
      lines(r_grad$Column,r_grad$beta_max,col="gray40")
      lines(r_grad$Column,r_grad$beta_median,col="gray40",lty=2)
      axis(1,at=r_grad$Column,cex.axis=1.2)
      if(is.na(yaxl)) axis(2,at=y_ax,cex.axis=1.2*(1-.2*(abs(y_ax)[1]>1000))) else axis(2,at=y_ax,labels = yaxl,cex.axis=1.2*(1-.2*(abs(y_ax)[1]>1000)))
      arrows(r_grad$Column, r_grad$beta_min, r_grad$Column, r_grad$beta_max+r_grad$sterr_max*1.96, length=0.05, angle=90, code=3,col="gray40")
      arrows(r_grad$Column, r_grad$beta_max, r_grad$Column, r_grad$beta_min-r_grad$sterr_min*1.96, length=0.05, angle=90, code=3,col="gray40")
    dev.off()
    
    tiff(paste0("Figures/Prop209_",v,"_Robustness_R2.tiff"),5,3.5,"in",res=300) #,width = 1000,height=1000
      plot(r_grad$Column,r_grad$r_min,
           ylim=c(0,max(r_grad$r_max)+0.01),
           pch=19, xlab="Specification: Number of Covariates", ylab="R-squared",axes=F,cex.lab=1.25,mgp=c(2.5,1,0))
      lines(r_grad$Column,r_grad$r_min,col="gray40")
      points(r_grad$Column,r_grad$r_max,pch=19)
      lines(r_grad$Column,r_grad$r_max,col="gray40")
      lines(r_grad$Column,r_grad$r_median,col="gray40",lty=2)
      axis(1,at=r_grad$Column,cex.axis=1.2)
      axis(2,at = seq(0,max(r_grad$r_max)+0.02,0.05),cex.axis=1.2)
    dev.off()
    print(paste(v,"  ",paste(r_grad$r_median[c(1,2,nrow(r_grad))],collapse="   ")))
  }
  
  
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces Figures IV, VI, A-4, A-8, A-11, and A-16.
  ##########################################################################################################
  ##########################################################################################################
  #Annual diff-in-diffs
  for(t in c("All_Fr","Asian")){
    A$urm<-A$URM ; if(t=="Asian"){
      A$urm[!A$CETHNICA%in%c("P")]<-NA
      A$urm[A$CETHNICA%in%c("D","E","G","H","L","N","V")]<-T
    }
    for(y in 1994:2002){
      A[,paste0("URM_",y)]<-A$urm*(A$YEARAPAY==y)
      A[,paste0("NonURM_",y)]<-(!A$urm)*(A$YEARAPAY==y)
      A[,paste0("Black_",y)]<-A$CETHNICA%in%c("B")*(A$YEARAPAY==y)*ifelse(A$CETHNICA%in%"A",NA,1)
      A[,paste0("Hispanic_",y)]<-A$CETHNICA%in%c("C","J")*(A$YEARAPAY==y)
    }
    
    
    d<-data.frame(t=c(1994:1996,1998:2002,1997),t95=c(1995:1996,1998:2002,1997,NA))
    d_single<-data.frame(t=c(1994:1996,1998:2002,1997),t95=c(1995:1996,1998:2002,1997,NA))
    d_multiurm<-data.frame(t=c(1994:1996,1998:2002,1997),t95=c(1995:1996,1998:2002,1997,NA))
    for(i in c(1,3:9)){
      A$DepVar<-A[,paste0("ADM0",i)]
      hi<-summary(felm(DepVar~URM_1994+URM_1995+URM_1996+URM_1998+URM_1999+URM_2000+URM_2001+URM_2002+factor(CETHNICA) + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY,data=A[A$YEARAPAY%in%1994:2002,]),robust=T)$coefficients[1:8,1:2]
      d[,paste0("c",i)]<-c(hi[,1],0)*100 ; d[,paste0("sd",i)]<-c(hi[,2],0.0000001)*100
      
      A$DepVar<-A[,paste0("ENR0",i)]
      hi<-summary(felm(DepVar~URM_1994+URM_1995+URM_1996+URM_1998+URM_1999+URM_2000+URM_2001+URM_2002+factor(CETHNICA) + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY,data=A[A$YEARAPAY%in%1994:2002,]),robust=T)$coefficients[1:8,1:2]
      d[,paste0("Ec",i)]<-c(hi[,1],0)*100 ; d[,paste0("Esd",i)]<-c(hi[,2],0.0000001)*100
    }
    for(i in c("Total_Wages_logCondit","Total_Wages_Rank_Condit","Total_Wages_Rank_Fixed_Condit","Total_Wages30s_Num100000")){
      mult<-100 ; if(grepl("Total_Wages",i)) mult<-1
      A$DepVar<-A[,i]
      hi<-summary(felm(DepVar~URM_1994+URM_1995+URM_1996+URM_1998+URM_1999+URM_2000+URM_2001+URM_2002+factor(CETHNICA) + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY,data=A[A$YEARAPAY%in%1994:2002,]),robust=T)$coefficients[1:8,1:2]
      d[,paste0("c",i)]<-c(hi[,1],0)*mult ; d[,paste0("sd",i)]<-c(hi[,2],0.0000001)*mult
      
      #Single-diff
      hi<-summary(felm(DepVar~URM_1994+URM_1995+URM_1996+URM_1998+URM_1999+URM_2000+URM_2001+URM_2002+NonURM_1994+NonURM_1995+NonURM_1996+NonURM_1998+NonURM_1999+NonURM_2000+NonURM_2001+NonURM_2002+factor(CETHNICA) + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom),data=A[A$YEARAPAY%in%1994:2002,]),robust=T)$coefficients[1:16,1:2]
      d_single[,paste0("c",i)]<-c(hi[1:8,1],0)*mult ; d_single[,paste0("sd",i)]<-c(hi[1:8,2],0.0000001)*mult ; d_single[,paste0("c",i,"_NonURM")]<-c(hi[9:16,1],0)*mult ; d_single[,paste0("sd",i,"_NonURM")]<-c(hi[9:16,2],0.0000001)*mult
      
      #Multi_URM
      hi<-summary(felm(DepVar~Black_1994+Black_1995+Black_1996+Black_1998+Black_1999+Black_2000+Black_2001+Black_2002+Hispanic_1994+Hispanic_1995+Hispanic_1996+Hispanic_1998+Hispanic_1999+Hispanic_2000+Hispanic_2001+Hispanic_2002+factor(CETHNICA) + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY,data=A[A$YEARAPAY%in%1994:2002,]),robust=T)$coefficients[1:16,1:2]
      d_multiurm[,paste0("c",i,"_Black")]<-c(hi[1:8,1],0)*mult ; d_multiurm[,paste0("sd",i,"_Black")]<-c(hi[1:8,2],0.0000001)*mult
      d_multiurm[,paste0("c",i,"_Hispanic")]<-c(hi[9:16,1],0)*mult ; d_multiurm[,paste0("sd",i,"_Hispanic")]<-c(hi[9:16,2],0.0000001)*mult
    }
    for(i in c("FE_WageVA_Chetty_AA","UnivGrad_SixYears","STEM_UCNSC","MA")){ #Drop 1994 for NSC
      mult<-100 ; if(grepl("FE",i)) mult<-1
      A$DepVar<-A[,i]
      if(i%in%c("UnivGrad_SixYears")){ ; A$DepVar[NA_to_T(A$S_Q!=1)]<-NA ; }
      hi<-summary(felm(DepVar~URM_1995+URM_1996+URM_1998+URM_1999+URM_2000+URM_2001+URM_2002+factor(CETHNICA) + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY,data=A[A$YEARAPAY%in%1995:2002,]),robust=T)$coefficients[1:7,1:2]
      d[,paste0("c",i)]<-c(hi[,1],0,NA)*mult ; d[,paste0("sd",i)]<-c(hi[,2],0.0000001,NA)*mult
      
      #Single-diff
      hi<-summary(felm(DepVar~URM_1995+URM_1996+URM_1998+URM_1999+URM_2000+URM_2001+URM_2002+NonURM_1995+NonURM_1996+NonURM_1998+NonURM_1999+NonURM_2000+NonURM_2001+NonURM_2002+factor(CETHNICA) + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom),data=A[A$YEARAPAY%in%1995:2002,]),robust=T)$coefficients[1:14,1:2]
      d_single[,paste0("c",i)]<-c(hi[1:7,1],0,NA)*mult ; d_single[,paste0("sd",i)]<-c(hi[1:7,2],0.0000001,NA)*mult ; d_single[,paste0("c",i,"_NonURM")]<-c(hi[8:14,1],0,NA)*mult ; d_single[,paste0("sd",i,"_NonURM")]<-c(hi[8:14,2],0.0000001,NA)*mult
      
      #Multi-URM
      hi<-summary(felm(DepVar~Black_1995+Black_1996+Black_1998+Black_1999+Black_2000+Black_2001+Black_2002+Hispanic_1995+Hispanic_1996+Hispanic_1998+Hispanic_1999+Hispanic_2000+Hispanic_2001+Hispanic_2002+factor(CETHNICA) + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY,data=A[A$YEARAPAY%in%1995:2002,]),robust=T)$coefficients[1:14,1:2]
      d_multiurm[,paste0("c",i,"_Black")]<-c(hi[1:7,1],0,NA)*mult ; d_multiurm[,paste0("sd",i,"_Black")]<-c(hi[1:7,2],0.0000001,NA)*mult 
      d_multiurm[,paste0("c",i,"_Hispanic")]<-c(hi[8:14,1],0,NA)*mult ; d_multiurm[,paste0("sd",i,"_Hispanic")]<-c(hi[8:14,2],0.0000001,NA)*mult
    }
    d<-d[order(d$t),]
    
    for(k in c("","E")){
      name<-"Adm" ; if(k=="E") name<-"Enr"
      lname<-"Admission" ; if(k=="E") lname<-"Enrollment"
      ylim<-c(-35,20) ; ylab<-c("-30","-15","0","15") ; yloc<-seq(-30,15,15)
      if(k=="E"){
        ylim<-c(-7,6) ; ylab<-c("-6","-3","0","3","6") ; yloc<-seq(-6,6,3)
      } 
      if(t=="Asian") ylim<-ylim/2
      
      tiff(paste0("Figures/AA_High",name,"_",t,".tiff"),5,3.5,"in",res=300)
        plot(d$t-0.15, d[,paste0(k,"c",1)],
             ylim=ylim,xlim=c(1993.6,2002.4),
             pch=16, xlab="", ylab=paste0("Percent"),cex.lab=1.1,cex.axis=1.2,mgp=c(2.4,1,0),yaxt="n")
        box()
        axis(2, labels = ylab, at = yloc,cex.axis=1.05)
        lines(c(1997.5,1997.5),c(-100,100),col="gray60")
        lines(c(1995.5,1995.5),c(-100,100),col="gray80",lty=2)
        lines(c(0,2100),rep(0,2),col="gray")
        i<-1 ; points(d$t-0.15, d[,paste0(k,"c",i)],pch=16,col="blue")
        arrows(d$t-.15, d[,paste0(k,"c",i)]-d[,paste0(k,"sd",i)]*1.96, d$t-.15, d[,paste0(k,"c",i)]+d[,paste0(k,"sd",i)]*1.96, length=0.05, angle=90, code=3,col="blue")
        i<-4 ; points(d$t, d[,paste0(k,"c",i)],pch=18,col="chartreuse4")
        arrows(d$t, d[,paste0(k,"c",i)]-d[,paste0(k,"sd",i)]*1.96, d$t, d[,paste0(k,"c",i)]+d[,paste0(k,"sd",i)]*1.96, length=0.05, angle=90, code=3,col="chartreuse4")
        i<-6 ; points(d$t+.15, d[,paste0(k,"c",i)],pch=15,col="firebrick1")
        arrows(d$t+.15, d[,paste0(k,"c",i)]-d[,paste0(k,"sd",i)]*1.96, d$t+.15, d[,paste0(k,"c",i)]+d[,paste0(k,"sd",i)]*1.96, length=0.05, angle=90, code=3,col="firebrick1")
        legend(x=1998, y=(23-16.7*(k=="E")-25*(t=="Asian")+19*(t=="Asian")*(k=="E")), c("Berkeley","UCLA","San Diego"), pch=c(16,18,15), col=c("blue","chartreuse4","firebrick1"),box.lty = 0,cex=.85,bty="n")
      dev.off()
      tiff(paste0("Figures/AA_Mid",name,"_",t,".tiff"),5,3.5,"in",res=300)
        plot(d$t-.15, d[,paste0(k,"c",3)],
             ylim=ylim,xlim=c(1993.6,2002.4),
             pch=16, xlab="", ylab=paste0("Percent"),cex.lab=1.1,cex.axis=1.2,mgp=c(2.4,1,0),yaxt="n")
        box()
        axis(2, labels = ylab, at = yloc,cex.axis=1.05)
        lines(c(1997.5,1997.5),c(-100,100),col="gray60")
        lines(c(1995.5,1995.5),c(-100,100),col="gray80",lty=2)
        lines(c(0,2100),rep(0,2),col="gray")
        i<-3 ; points(d$t-.15, d[,paste0(k,"c",i)],pch=16,col="blue")
        arrows(d$t-.15, d[,paste0(k,"c",i)]-d[,paste0(k,"sd",i)]*1.96, d$t-.15, d[,paste0(k,"c",i)]+d[,paste0(k,"sd",i)]*1.96, length=0.05, angle=90, code=3,col="blue")
        i<-8 ; points(d$t, d[,paste0(k,"c",i)],pch=18,col="chartreuse4")
        arrows(d$t, d[,paste0(k,"c",i)]-d[,paste0(k,"sd",i)]*1.96, d$t, d[,paste0(k,"c",i)]+d[,paste0(k,"sd",i)]*1.96, length=0.05, angle=90, code=3,col="chartreuse4")
        i<-9 ; points(d$t+.15, d[,paste0(k,"c",i)],pch=15,col="firebrick1") #Color 4: darkorange
        arrows(d$t+.15, d[,paste0(k,"c",i)]-d[,paste0(k,"sd",i)]*1.96, d$t+.15, d[,paste0(k,"c",i)]+d[,paste0(k,"sd",i)]*1.96, length=0.05, angle=90, code=3,col="firebrick1")
        legend(x=1998, y=23.5-25*(k=="E")-27*(t=="Asian")+27.25*(t=="Asian")*(k=="E"), c("Davis","S. Barbara","Irvine"), pch=c(16,18,15), col=c("blue","chartreuse4","firebrick1"),box.lty = 0,cex=.85,bty="n")
      dev.off()
      tiff(paste0("Figures/AA_Low",name,"_",t,".tiff"),5,3.5,"in",res=300)
        plot(d$t-.15, d[,paste0(k,"c",7)],
             ylim=ylim,xlim=c(1993.6,2002.4),
             pch=16, xlab="", ylab=paste0("Percent"),cex.lab=1.1,cex.axis=1.2,mgp=c(2.4,1,0),yaxt="n")
        box()
        axis(2, labels = ylab, at = yloc,cex.axis=1.05)
        lines(c(1997.5,1997.5),c(-100,100),col="gray60")
        lines(c(1995.5,1995.5),c(-100,100),col="gray80",lty=2)
        lines(c(0,2100),rep(0,2),col="gray")
        i<-7 ; points(d$t-.15, d[,paste0(k,"c",i)],pch=16,col="blue")
        arrows(d$t-.15, d[,paste0(k,"c",i)]-d[,paste0(k,"sd",i)]*1.96, d$t-.15, d[,paste0(k,"c",i)]+d[,paste0(k,"sd",i)]*1.96, length=0.05, angle=90, code=3,col="blue")
        i<-5 ; points(d$t, d[,paste0(k,"c",i)],pch=18,col="chartreuse4")
        arrows(d$t, d[,paste0(k,"c",i)]-d[,paste0(k,"sd",i)]*1.96, d$t, d[,paste0(k,"c",i)]+d[,paste0(k,"sd",i)]*1.96, length=0.05, angle=90, code=3,col="chartreuse4")
        legend(x=1998, y=20-21*(k=="E")-29*(t=="Asian")+29.25*(t=="Asian")*(k=="E"), c("Santa Cruz","Riverside"), pch=c(16,18), col=c("blue","chartreuse4"),box.lty = 0,cex=.85)
      dev.off()
    }
    
    o<-c("FE_WageVA_Chetty_AA","UnivGrad_SixYears","STEM_UCNSC","MA","Total_Wages_logCondit","Total_Wages_Rank_Condit","Total_Wages_Rank_Fixed_Condit","Total_Wages30s_Num100000")
    oname<-c("Dollars ('000s)","Percent","Percent","Percent","Log Dollars","Percentiles","Percentiles","Percent")
    for(i in o){
      if(i%in%c("FE_WageVA_Chetty_AA","UnivGrad_SixYears","STEM_UCNSC","MA")){ ; d$tt<-d$t95 ; d_single$tt<-d_single$t95 ; d_multiurm$tt<-d_multiurm$t95 } else{ ; d$tt<-d$t ; d_single$tt<-d_single$t ; d_multiurm$tt<-d_multiurm$t } #Use 1995 starting point for NSC
      d<-d[order(d$tt),] ; d_single<-d_single[order(d_single$tt),] ; d_multiurm<-d_multiurm[order(d_multiurm$tt),]
      ylim<-c(-10,10) ; yax<-seq(-10,10,10) ; yaxl<-seq(-10,10,10)
      if(t=="All_Fr"&i%in%c("UnivGrad_SixYears")){ ; ylim<-c(-10,3) ; yax<-seq(-10,2.5,2.5) ; yaxl<-c("-10",NA,"-5",NA,0,NA) ; }
      if(t=="All_Fr"&i%in%c("MA")){ ; ylim<-c(-4,2.5) ; yax<-seq(-4,2,2) ; yaxl<-c("-4","-2","0","2") ; }
      if(t=="All_Fr"&i%in%c("STEM_UCNSC")){ ; ylim<-c(-4,2.5) ; yax<-seq(-4,2,2) ; yaxl<-c("-4","-2","0","2") ; }
      if(t=="All_Fr"&i%in%c("Total_Wages_Rank_Condit")){ ; ylim<-c(-3,1.5) ; yax<-seq(-3,1) ; yaxl<-c("-3","-2","-1","0","1") ; }
      if(t=="All_Fr"&i%in%c("Total_Wages_Rank_Fixed_Condit")){ ; ylim<-c(-4.5,2.5) ; yax<-seq(-4,2,2) ; yaxl<-c("-4","-2","0","1") ; }
      if(t=="All_Fr"&i%in%c("Total_Wages30s_Num100000")){ ; ylim<-c(-.22,.2) ; yax<-seq(-.2,.2,.1) ; yaxl<-c("-0.2",NA,0,NA,0.2) ; }
      if(i%in%c("Total_Wages_logCondit")){ ; ylim<-c(-.15,.1) ; yax<-seq(-.15,.1,.05) ; yaxl<-c(NA,"-0.1",NA,0,NA,"0.1") ; }
      if(i%in%c("FE_WageVA_Chetty_AA")){ ; ylim<-c(-2000,1000) ; yax<-seq(-2000,1000,1000) ; yaxl<-c("-2","-1","0","1") ; }
      tiff(paste0("Figures/AA_",i,"_",t,".tiff"),5,3.5,"in",res=300)
        plot(d$tt, d[,paste0('c',i)],
             ylim=ylim,
             pch=16, xlab="", ylab=oname[o==i],cex.lab=1.1,cex.axis=1.2,mgp=c(2.4,1,0),yaxt="n"
        )
        axis(2, labels =yaxl , at = yax,cex.axis=1.1)
        lines(c(1997.5,1997.5),c(-100000,100000),col="gray60")
        lines(c(1995.5,1995.5),c(-100000,100000),col="gray80",lty=2)
        lines(c(0,2100),rep(0,2),col="gray")
        points(d$tt, d[,paste0("c",i)],pch=16)
        arrows(d$tt, d[,paste0("c",i)]-d[,paste0("sd",i)]*1.96, d$tt, d[,paste0("c",i)]+d[,paste0("sd",i)]*1.96, length=0.05, angle=90, code=3)
      dev.off()
      ylim<-ylim*2
      if(t=="All_Fr"&i%in%c("UnivGrad_SixYears")) ylim<-c(-10,10)
      if(t=="All_Fr"&i%in%c("MA")) ylim<-c(-8,2)
      if(t=="All_Fr"&i%in%c("Total_Wages_logCondit")) ylim<-c(-.2,.12)
      if(t=="All_Fr"&i%in%c("Total_Wages30s_Num100000")) ylim<-c(-.3,.2)
      if(i%in%c("FE_WageVA_Chetty_AA")) ylim<-c(-2500,1250)
      if(i%in%c("Total_Wages_Rank_Condit")) ylim<-c(-10,10)
      tiff(paste0("Figures/AA_",i,"_",t,"_SingleDiff.tiff"),5,3.5,"in",res=300)
        plot(d_single$tt-.1, d_single[,paste0('c',i)],
             ylim=ylim,xlim=c(min(d$tt,na.rm=T)-.1,max(d$tt,na.rm=T)+.1),
             pch=16, xlab="", ylab=oname[o==i],cex.lab=1.25,cex.axis=1.5
        )
        lines(c(1997.5,1997.5),c(-10000,10000),col="gray60")
        lines(c(1995.5,1995.5),c(-10000,10000),col="gray80",lty=2)
        lines(c(0,2100),rep(0,2),col="gray")
        points(d_single$tt-.1, d_single[,paste0("c",i)],pch=16)
        arrows(d_single$tt-.1, d_single[,paste0("c",i)]-d_single[,paste0("sd",i)]*1.96, d_single$tt-.1, d_single[,paste0("c",i)]+d_single[,paste0("sd",i)]*1.96, length=0.05, angle=90, code=3)
        points(d_single$tt+.1, d_single[,paste0("c",i,"_NonURM")],pch=16,col="gray50")
        arrows(d_single$tt+.1, d_single[,paste0("c",i,"_NonURM")]-d_single[,paste0("sd",i,"_NonURM")]*1.96, d_single$tt+.1, d_single[,paste0("c",i,"_NonURM")]+d_single[,paste0("sd",i,"_NonURM")]*1.96, length=0.05, angle=90, code=3,col="gray50")
      dev.off()
      if(t=="All_Fr"){
        if(i%in%c("FE_WageVA_Chetty_AA")){ ; ylim<-c(-2500,1000) ; yax<-seq(-2000,1000,1000) ; yaxl<-c("-2","-1","0","1") ; }
        if(t=="All_Fr"&i%in%c("UnivGrad_SixYears")){ ; ylim<-c(-12.5,4) ; yax<-seq(-10,2.5,2.5) ; yaxl<-c("-10",NA,"-5",NA,0,NA) ; }
        if(t=="All_Fr"&i%in%c("STEM_UCNSC")){ ; ylim<-c(-6,2) ; yax<-seq(-6,2,2) ; yaxl<-c("-6","-4","-2","0","2") ; }
        if(t=="All_Fr"&i%in%c("MA")){ ; ylim<-c(-8,4) ; yax<-seq(-8,4,4) ; yaxl<-c("-8","-4","0","4") ; }
        if(i%in%c("Total_Wages_logCondit")){ ; ylim<-c(-.15,.15) ; yax<-seq(-.15,.15,.05) ; yaxl<-c(NA,"-0.1",NA,0,NA,"0.1",NA) ; }
        if(t=="All_Fr"&i%in%c("Total_Wages_Rank_Condit")){ ; ylim<-c(-4,2.5) ; yax<-seq(-4,2,2) ; yaxl<-c("-4","-2","0","2") ; }
        if(t=="All_Fr"&i%in%c("Total_Wages_Rank_Fixed_Condit")){ ; ylim<-c(-5,3) ; yax<-seq(-4,2,2) ; yaxl<-c("-4","-2","0","2") ; }
        tiff(paste0("Figures/AA_",i,"_",t,"_MultiURM.tiff"),5,3.5,"in",res=300)
          plot(d_multiurm$tt-.1, d_multiurm[,paste0('c',i,'_Black')],
               ylim=ylim,xlim=c(min(d$tt,na.rm=T)-.1,max(d$tt,na.rm=T)+.1),
               pch=17, xlab="", ylab=oname[o==i],cex.lab=1.1,cex.axis=1.2,mgp=c(2.4,1,0),yaxt="n"
          )
          axis(2, labels =yaxl , at = yax,cex.axis=1.1)
          lines(c(1997.5,1997.5),c(-10000,10000),col="gray60")
          lines(c(1995.5,1995.5),c(-10000,10000),col="gray80",lty=2)
          lines(c(0,2100),rep(0,2),col="gray")
          points(d_multiurm$tt-.1, d_multiurm[,paste0("c",i,"_Black")],pch=17)
          arrows(d_multiurm$tt-.1, d_multiurm[,paste0("c",i,"_Black")]-d_multiurm[,paste0("sd",i,"_Black")]*1.96, d_multiurm$tt-.1, d_multiurm[,paste0("c",i,"_Black")]+d_multiurm[,paste0("sd",i,"_Black")]*1.96, length=0.05, angle=90, code=3)
          points(d_multiurm$tt+.1, d_multiurm[,paste0("c",i,"_Hispanic")],pch=15,col="gray50")
          arrows(d_multiurm$tt+.1, d_multiurm[,paste0("c",i,"_Hispanic")]-d_multiurm[,paste0("sd",i,"_Hispanic")]*1.96, d_multiurm$tt+.1, d_multiurm[,paste0("c",i,"_Hispanic")]+d_multiurm[,paste0("sd",i,"_Hispanic")]*1.96, length=0.05, angle=90, code=3,col="gray50")
        dev.off()
      }
    }
  }
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces Figure A-10.
  ##########################################################################################################
  ##########################################################################################################
  A$urm<-A$URM
  for(y in 1994:2002){
    A[,paste0("URM_",y)]<-A$urm*(A$YEARAPAY==y)
    A[,paste0("NonURM_",y)]<-(!A$urm)*(A$YEARAPAY==y)
    A[,paste0("Black_",y)]<-A$CETHNICA%in%c("B")*(A$YEARAPAY==y)*ifelse(A$CETHNICA%in%"A",NA,1)
    A[,paste0("Hispanic_",y)]<-A$CETHNICA%in%c("C","J")*(A$YEARAPAY==y)
  }
  reg_totalwages<-felm(Total_Wages_logCondit~URM_1994+URM_1995+URM_1996+URM_1998+URM_1999+URM_2000+URM_2001+URM_2002+factor(CETHNICA) + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY,data=A[A$YEARAPAY%in%1994:2002,])
  reg_rankwages<-felm(Total_Wages_Rank_Condit~URM_1994+URM_1995+URM_1996+URM_1998+URM_1999+URM_2000+URM_2001+URM_2002+factor(CETHNICA) + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY,data=A[A$YEARAPAY%in%1994:2002,])
  
  for(v in list(list("Total_Wages_logCondit",reg_totalwages,0.015,3),
                list("Total_Wages_Rank_Condit",reg_rankwages,0.3,1))){
    l_vec<-c(1,1,0,0,0)
    reg<-v[[2]]
    Mvec<-seq(from = 0, to = v[[3]], by = v[[3]]/8)
    hi<-createSensitivityResults(betahat=reg$coefficients[1:8],sigma=reg$vcv[1:8,1:8],numPrePeriods=3,numPostPeriods=5,l_vec=l_vec,Mvec=Mvec,parallel = T)
    orig<-constructOriginalCS(betahat=reg$coefficients[1:8],sigma=reg$vcv[1:8,1:8],numPrePeriods=3,numPostPeriods=5,l_vec=l_vec)
    scaleFUN <- function(x) sprintf(paste0("%.",v[[4]],"f"), x)
    g<-createSensitivityPlot(hi,orig)+
      theme(axis.text.y   = element_text(size=17.5),
            axis.text.x   = element_text(size=17.5),
            axis.title.y  = element_text(size=17.5),
            axis.title.x  = element_text(size=17.5),
            panel.background = element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),
            panel.border = element_rect(colour = "black", fill=NA, size=1)
      )+ scale_x_continuous(labels=scaleFUN)
    ggsave(filename=paste0("RRSens_",v[[1]],".tiff"),g,"tiff",path = "Figures",width = 5,height=3.5,units = "in",dpi = 300)
  }
  
  
  #Show non-parametric outcome changes by AI
  A$UCBEnrFirstIncCC<-NA_to_F(A$NSC_EnrFirstIncCC_College_Code%in%c(1312))&NA_to_F(!grepl("UNIVERSITY OF CALIFORNIA-EXTENSION",A$NSC_EnrFirstIncCC_College_Name)) 
  A$UCLAEnrFirstIncCC<-NA_to_F(A$NSC_EnrFirstIncCC_College_Code%in%c(1315))&NA_to_F(!grepl("UNIVERSITY OF CALIFORNIA-EXTENSION",A$NSC_EnrFirstIncCC_College_Name)) 
  A$UCSDEnrFirstIncCC<-NA_to_F(A$NSC_EnrFirstIncCC_College_Code%in%c(1317))&NA_to_F(!grepl("UNIVERSITY OF CALIFORNIA-EXTENSION",A$NSC_EnrFirstIncCC_College_Name)) 
  A$UCSBEnrFirstIncCC<-NA_to_F(A$NSC_EnrFirstIncCC_College_Code%in%c(1320))&NA_to_F(!grepl("UNIVERSITY OF CALIFORNIA-EXTENSION",A$NSC_EnrFirstIncCC_College_Name))
  A$UCIEnrFirstIncCC<-NA_to_F(A$NSC_EnrFirstIncCC_College_Code%in%c(1314))&NA_to_F(!grepl("UNIVERSITY OF CALIFORNIA-EXTENSION",A$NSC_EnrFirstIncCC_College_Name)) 
  A$UCDEnrFirstIncCC<-NA_to_F(A$NSC_EnrFirstIncCC_College_Code%in%c(1313))&NA_to_F(!grepl("UNIVERSITY OF CALIFORNIA-EXTENSION",A$NSC_EnrFirstIncCC_College_Name)) 
  A$UCSCEnrFirstIncCC<-NA_to_F(A$NSC_EnrFirstIncCC_College_Code%in%c(1321))&NA_to_F(!grepl("UNIVERSITY OF CALIFORNIA-EXTENSION",A$NSC_EnrFirstIncCC_College_Name)) 
  A$UCREnrFirstIncCC<-NA_to_F(A$NSC_EnrFirstIncCC_College_Code%in%c(1316))&NA_to_F(!grepl("UNIVERSITY OF CALIFORNIA-EXTENSION",A$NSC_EnrFirstIncCC_College_Name)) 
  A$UCBUCLAEnrFirstIncCC<-NA_to_F(A$NSC_EnrFirstIncCC_College_Code%in%c(1312,1315))&NA_to_F(!grepl("UNIVERSITY OF CALIFORNIA-EXTENSION",A$NSC_EnrFirstIncCC_College_Name))
  A$MidUCSDEnrFirstIncCC<-NA_to_F(A$NSC_EnrFirstIncCC_College_Code%in%c(1317,1320,1314,1313))&NA_to_F(!grepl("UNIVERSITY OF CALIFORNIA-EXTENSION",A$NSC_EnrFirstIncCC_College_Name))
  A$CAPrivIncStanEnrFirstIncCC<-NA_to_F(A$NSC_EnrFirstIncCC_College_Code%in%c(1305))|A$CAPrivEnrFirstIncCC
  A$NonCAIncIvyEnrFirstIncCC<-!NA_to_F(A$NSC_EnrFirstIncCC_College_Code%in%c(1305))&(A$NonCAEnrFirstIncCC|A$IvyEnrFirstIncCC)
  A$AppMid<-(A$APP03==1|A$APP06==1|A$APP08==1|A$APP09==1)
  A$AppLow<-(A$APP07==1|A$APP05==1)
  A$ADMMid<-(A$ADM03==1|A$ADM06==1|A$ADM08==1|A$ADM09==1) ; A$ADMMin[A$AppMid==0]<-NA
  A$ADMLow<-(A$ADM07==1|A$ADM05==1) ; A$ADMLow[A$AppLow==0]<-NA
  A$ADMVeryHigh<-(A$ADM01==1|A$ADM04==1) ; A$ADMHigh[!(A$APP01==1|A$APP04==1)]<-NA
  A$UCEnrFirstIncCC<-(A$HighUCEnrFirstIncCC==1|A$MidUCEnrFirstIncCC==1|A$LowUCEnrFirstIncCC==1)
  vars<-c("UCBEnrFirstIncCC","UCLAEnrFirstIncCC","UCSDEnrFirstIncCC","UCSBEnrFirstIncCC","UCIEnrFirstIncCC","UCDEnrFirstIncCC","UCSCEnrFirstIncCC","UCREnrFirstIncCC","HighUCEnrFirstIncCC","MidUCEnrFirstIncCC","LowUCEnrFirstIncCC","UCEnrFirstIncCC","CSUEnrFirstIncCC","CCEnrFirstIncCC","IvyEnrFirstIncCC","CAPrivEnrFirstIncCC","NonCAEnrFirstIncCC","NonEnrFirstIncCC","UCBUCLAEnrFirstIncCC","MidUCSDEnrFirstIncCC","CAPrivIncStanEnrFirstIncCC","NonCAIncIvyEnrFirstIncCC",
          "FE_GradVA_Chetty_AA","FE_WageVA_Chetty_AA",
          "UnivGrad_SixYears","MA","MA_STEM",
          "Total_Wages_Condit","Total_Wages30s_Condit","Total_Wages_logCondit","Total_Wages30s_Num100000","wage_sum_15",
          "APP01","APP03","APP04","APP05","APP06","APP07","APP08","APP09","AppMid","AppLow",
          "ADM01","ADM03","ADM04","ADM05","ADM06","ADM07","ADM08","ADM09","ADMMid","ADMLow","ADMVeryHigh")
  A$Post<-A$YEARAPAY>1997
  
  b<-unique(c(quantile(A$S[A$YEARAPAY%in%1996:1999&A$URM], probs = seq(0, 1, by = 1/100),na.rm=T)))
  A$Bin<-as.integer(cut(A$S, breaks=b,labels=1:(length(b)-1), include.lowest=TRUE))
  x1<-A[A$YEARAPAY%in%1996:1999,]
  x1$FE_WageVA_Chetty_AA[is.na(x1$Total_Wages30s_Condit)]<-NA
  hi_bin<-data.frame()
  for(i in 3:98){
    x<-x1[x1$Bin%in%(i-14):(i+14),]
    x$weight<- 1-(abs(x$Bin-i)/15)
    y<-data.frame(Post=c(0,0,1,1),URM=c(0,1,0,1),Bin=rep(i,4))
    for(v in c(vars)) y[,v]<-c(weighted.mean(x[x$Post==0&x$URM==0,v],x$weight[x$Post==0&x$URM==0],na.rm=T),
                                                     weighted.mean(x[x$Post==0&x$URM==1,v],x$weight[x$Post==0&x$URM==1],na.rm=T),
                                                     weighted.mean(x[x$Post==1&x$URM==0,v],x$weight[x$Post==1&x$URM==0],na.rm=T),
                                                     weighted.mean(x[x$Post==1&x$URM==1,v],x$weight[x$Post==1&x$URM==1],na.rm=T))
    hi_bin<-rbind(hi_bin,y)
    print(i)
  }
  hi_bin$UnivGrad_SixYears<-hi_bin$UnivGrad_SixYears*100
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces Figures II, III, and A-7.
  ##########################################################################################################
  ##########################################################################################################
  for(v in c(vars)){
    if(sum(inrange(hi_bin[,v],-1,1))==nrow(hi_bin)){
      if(0==1) hi[,v]<-hi[,v]*100
      hi_bin[,v]<-hi_bin[,v]*100
    }
    ylab<-"Percentage Points" ; if(grepl("Wages",v)&!grepl("Num",v)){ ; if(grepl("log",v)) ylab<-"Log Dollars" else ylab<-"Dollars" ; }
    ylim<-c(min(hi_bin[hi_bin$Post&hi_bin$URM,v]-hi_bin[!hi_bin$Post&hi_bin$URM,v]),max(hi_bin[hi_bin$Post&hi_bin$URM,v]-hi_bin[!hi_bin$Post&hi_bin$URM,v]))
    if(grepl("Enr",v)) ylim<-c(-10,8)
    if(grepl("A[Pp][Pp]",v)) ylim<-c(-7,15)
    if(grepl("ADM",v)) ylim<-c(-65,5)
    ylimdd<-c(min(hi_bin[hi_bin$Post&hi_bin$URM,v]-hi_bin[!hi_bin$Post&hi_bin$URM,v]-(hi_bin[hi_bin$Post&!hi_bin$URM,v]-hi_bin[!hi_bin$Post&!hi_bin$URM,v])),max(hi_bin[hi_bin$Post&hi_bin$URM,v]-hi_bin[!hi_bin$Post&hi_bin$URM,v]-(hi_bin[hi_bin$Post&!hi_bin$URM,v]-hi_bin[!hi_bin$Post&!hi_bin$URM,v])))
    if(grepl("Enr",v)) ylimdd<-c(-16,16)
    if(grepl("A[Pp][Pp]",v)) ylimdd<-c(-7,15)
    if(grepl("ADM",v)) ylimdd<-c(-65,5)
    ylimnod<-c(0,50)
    if(v=="ADMVeryHigh") ylimnod<-c(0,100)
    if(v%in%c("UCBUCLAEnrFirstIncCC","UCEnrFirstIncCC")) ylimnod<-c(0,70)
    
    #Quantile
    tiff(paste0("Figures/ResPlot_NoD_",v,"_AIQuant.tiff"),5,3.5,"in",res=300)
      plot(hi_bin$Bin[hi_bin$Post&hi_bin$URM],hi_bin[!hi_bin$Post&hi_bin$URM,v],type="l",mgp=c(2.5,1,0),pch=NA,xlab="AI Percentile",ylab=ylab,ylim=ylimnod,cex.lab=1.25,cex.axis=1.2,lwd=2) #Single-Diff
      lines(hi_bin$Bin[hi_bin$Post&hi_bin$URM],hi_bin[hi_bin$Post&hi_bin$URM,v],lty=2,col="black",lwd=2)
      lines(hi_bin$Bin[hi_bin$Post&hi_bin$URM],hi_bin[!hi_bin$Post&!hi_bin$URM,v],lty=1,col="gray50",lwd=2)
      lines(hi_bin$Bin[hi_bin$Post&hi_bin$URM],hi_bin[hi_bin$Post&!hi_bin$URM,v],lty=2,col="gray50",lwd=2)
      lines(c(-10000,10000),c(0,0))  
      abline(v=25,col="gray50",lty=3)
      abline(v=50,col="gray50",lty=3)
      abline(v=75,col="gray50",lty=3)
    dev.off()
    tiff(paste0("Figures/ResPlot_SD_",v,"_AIQuant.tiff"),5,3.5,"in",res=300)
      plot(hi_bin$Bin[hi_bin$Post&hi_bin$URM],hi_bin[hi_bin$Post&hi_bin$URM,v]-hi_bin[!hi_bin$Post&hi_bin$URM,v],type="l",mgp=c(2.5,1,0),pch=NA,xlab="AI Percentile",ylab=ylab,ylim=ylim,cex.lab=1.25,cex.axis=1.2,lwd=2) #Single-Diff
      lines(hi_bin$Bin[hi_bin$Post&hi_bin$URM],hi_bin[hi_bin$Post&!hi_bin$URM,v]-hi_bin[!hi_bin$Post&!hi_bin$URM,v],lty=2,col="gray50",lwd=2)
      lines(c(-10000,10000),c(0,0))  
      abline(v=25,col="gray50",lty=3)
      abline(v=50,col="gray50",lty=3)
      abline(v=75,col="gray50",lty=3)
      if(grepl("UCB|UCLA",v)) legend(100,8.1,c("URM","Non-URM"),lty=c(1,2),col=c("black","gray50"),cex=.75,xjust=1,lwd=c(2,2)) else if(!grepl("A[pP][Pp]|ADM",v)) legend(100,-10.1,c("URM","Non-URM"),lty=c(1,2),col=c("black","gray50"),cex=.8,xjust=1,yjust=0,lwd=c(2,2))
    dev.off()
    tiff(paste0("Figures/ResPlot_DD_",v,"_AIQuant.tiff"),5,3.5,"in",res=300)
      plot(hi_bin$Bin[hi_bin$Post&hi_bin$URM],hi_bin[hi_bin$Post&hi_bin$URM,v]-hi_bin[!hi_bin$Post&hi_bin$URM,v]-(hi_bin[hi_bin$Post&!hi_bin$URM,v]-hi_bin[!hi_bin$Post&!hi_bin$URM,v]),type="l",pch=NA,xlab="AI Percentile",ylab=ylab,ylim=ylimdd,cex.lab=1.25,cex.axis=1.5,lwd=2) #Diff-in-Diff
      lines(c(-10000,10000),c(0,0),col="gray50")  
      abline(v=25,col="gray50",lty=2)
      abline(v=50,col="gray50",lty=2)
      abline(v=75,col="gray50",lty=2)
    dev.off()
  }
  tiff(paste0("Figures/Prop209_SD_Legend.tiff"),width=8,height=3,units="in",res=500)
    plot(100000,100000,bty="n",box=F,yaxt="n",xaxt="n",xlim=c(0,50),ylim=c(0,10),xlab="",ylab="")
    legend(1,8,legend=c("URM","Non-URM"),lty=c(1,2),lwd=c(1,1),col=c("black","gray50"),cex=.6,ncol = 2)
    dev.off()
    tiff(paste0("Figures/Prop209_NoD_Legend.tiff"),width=8,height=3,units="in",res=1000)
    plot(100000,100000,bty="n",box=F,yaxt="n",xaxt="n",xlim=c(0,100),ylim=c(0,10),xlab="",ylab="")
    legend(1,8,legend=c("URM Pre-209","URM Post-209","Non-URM Pre-209","Non-URM Post-209"),lty=c(1,2,1,2),lwd=c(1,1,1,1),col=c("black","black","gray50","gray50"),cex=.3,ncol = 4)
  dev.off()
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces Figure A-18.
  ##########################################################################################################
  ##########################################################################################################
  load("Figures/Tables.Rda")
  for(var in list(list("FE_GradVA_Chetty_AA","UnivGrad_SixYears","BA"),
                  list("FE_WageVA_Chetty_AA","Total_Wages30s_Condit","WageVA"))){
    v<-var[[1]]
    v1<-var[[2]]
    if(v1%in%c("UnivGrad_SixYears_NetHS","UnivGrad_SixYears")) ybar<-as.numeric(gsub("[^0-9.-]","",R2$UnivGrad_SixYears[4]))
    if(v1%in%c("Total_Wages30s_Condit_NetHS","Total_Wages30s_Condit")) ybar<-as.numeric(gsub("[^0-9.-]","",R2$Total_Wages30s_Condit[4]))
    ylab<-"Percentage Points" ; if(grepl("Wage",v)&!grepl("Num",v)){ ; if(grepl("log",v)) ylab<-"Log Dollars" else ylab<-"Dollars" ; }
    tiff(paste0("Figures/ResPlot_DDCombined_",var[[3]],"_AIQuant.tiff"),5,3.5,"in",res=300)
      plot(hi_bin$Bin[hi_bin$Post&hi_bin$URM],hi_bin[hi_bin$Post&hi_bin$URM,v1]-hi_bin[!hi_bin$Post&hi_bin$URM,v1]-(hi_bin[hi_bin$Post&!hi_bin$URM,v1]-hi_bin[!hi_bin$Post&!hi_bin$URM,v1]),type="l",pch=NA,xlab="Academic Index Percentile",ylab=ylab,cex.lab=1.15,cex.axis=1.2,lwd=2) #Diff-in-Diff
      lines(hi_bin$Bin[hi_bin$Post&hi_bin$URM],hi_bin[hi_bin$Post&hi_bin$URM,v]-hi_bin[!hi_bin$Post&hi_bin$URM,v]-(hi_bin[hi_bin$Post&!hi_bin$URM,v]-hi_bin[!hi_bin$Post&!hi_bin$URM,v]),col="blue",lty=2,lwd=2)
      lines(c(-10000,10000),c(0,0),col="gray50")  
      #abline(h=ybar,lty=3,col="gray60")   
      abline(v=25,col="gray50",lty=3)
      abline(v=50,col="gray50",lty=3)
      abline(v=75,col="gray50",lty=3)
    dev.off()
  }
  tiff(paste0("Figures/ResPlot_DDCombined_Legend.tiff"),width=8,height=3,units="in",res=500)
    plot(100000,100000,bty="n",box=F,yaxt="n",xaxt="n",xlim=c(0,50),ylim=c(0,10),xlab="",ylab="")
    legend(1,8,legend=c("CFSTY Value-Added","Observed Outcome"),lty=c(2,1),lwd=c(2,2),col=c("blue","black"),cex=.6,ncol = 2)
  dev.off()
  
  
  ###
  #  How many URM students exited CA?
  ###
  mean(A$IvyEnrFirstIncCC[A$URM==1&NA_to_F(A$S_Q==4)&A$YEARAPAY%in%1996:1997],na.rm=T)
  mean(A$IvyEnrFirstIncCC[A$URM==1&NA_to_F(A$S_Q==4)&A$YEARAPAY%in%1998:1999],na.rm=T)
  table(A$IvyEnrFirstIncCC[A$URM==1&A$S_Q==4],round(A$YEARAPAY/2-.4)[A$URM==1&A$S_Q==4])
  table(A$IvyEnrFirstIncCC[A$URM==1&A$S_Q==4],round(A$YEARAPAY/2-.4)[A$URM==1&A$S_Q==4])[2,]/table(round(A$YEARAPAY/2-.4)[A$URM==1&A$S_Q==4])
  table(A$NonCAIncIvyEnrFirstIncCC[A$URM==1&A$S_Q==2],A$YEARAPAY[A$URM==1&A$S_Q==2])[2,]/table(A$YEARAPAY[A$URM==1&A$S_Q==4])
  R<-data.frame(Time=6:16,Emp=NA,K100=NA)
  for(i in 6:16){
    R$Emp[i-5]<-mean(A[A$IvyEnrFirstIncCC==1&A$URM&A$YEARAPAY%in%1998:1999,paste0("wage_sum_",i,"_Employed")])*100
    R$EmpNon[i-5]<-mean(A[A$IvyEnrFirstIncCC==0&A$URM&A$YEARAPAY%in%1998:1999,paste0("wage_sum_",i,"_Employed")])*100
    R$EmpS4[i-5]<-mean(A[A$IvyEnrFirstIncCC==0&A$URM&A$YEARAPAY%in%1998:1999&NA_to_F(A$S_Q==4),paste0("wage_sum_",i,"_Employed")])*100
    R$Emp96[i-5]<-mean(A[A$IvyEnrFirstIncCC==1&A$URM&A$YEARAPAY%in%1996:1997,paste0("wage_sum_",i,"_Employed")])*100
    R$EmpNon96[i-5]<-mean(A[A$IvyEnrFirstIncCC==0&A$URM&A$YEARAPAY%in%1996:1997,paste0("wage_sum_",i,"_Employed")])*100
    R$EmpS496[i-5]<-mean(A[A$IvyEnrFirstIncCC==0&A$URM&A$YEARAPAY%in%1996:1997&NA_to_F(A$S_Q==4),paste0("wage_sum_",i,"_Employed")])*100
    R$K100[i-5]<-mean(A[A$IvyEnrFirstIncCC==1&A$URM&A$YEARAPAY%in%1998:1999,paste0("wage_sum_",i,"_100000")])*100
    R$K10096[i-5]<-mean(A[A$IvyEnrFirstIncCC==1&A$URM&A$YEARAPAY%in%1996:1997,paste0("wage_sum_",i,"_100000")])*100
    R$K100S4[i-5]<-mean(A[A$IvyEnrFirstIncCC==0&A$URM&A$YEARAPAY%in%1998:1999&NA_to_F(A$S_Q==4),paste0("wage_sum_",i,"_100000")])*100
    R$K100S496[i-5]<-mean(A[A$IvyEnrFirstIncCC==0&A$URM&A$YEARAPAY%in%1996:1997&NA_to_F(A$S_Q==4),paste0("wage_sum_",i,"_100000")])*100
  } 
  plot(R$Time,R$Emp,type="l",pch=NA,xlab="Years After Application",ylab="Proportion of Students",xlim=c(6,16),ylim=c(0,80),cex.lab=1.25,cex.axis=1.2,mgp=c(2.6,1,0),col="blue",lty=2) #Single-Diff
    lines(c(-10000,10000000),c(0,0))
    lines(R$Time,R$Emp96,col="blue")
    lines(R$Time,R$EmpNon,lty=2)
    lines(R$Time,R$EmpNon96)
    lines(R$Time,R$EmpS4,col="gray50",lty=2)
    lines(R$Time,R$EmpS496,col="gray50")
    lines(R$Time,R$K100/R$Emp*100,col="forestgreen",lwd=1.5,lty=2)
    lines(R$Time,R$K10096/R$Emp96*100,col="forestgreen",lwd=1.5)
    lines(R$Time,R$K100S4/R$EmpS4*100,col="red",lwd=1.5,lty=2)
    lines(R$Time,R$K100S496/R$EmpS496*100,col="red",lwd=1.5)
    text(x=12,y=76,pos=4,labels="Non-Ivy+ CA Emp.",cex = .6,font=2)
    text(x=5.5,y=66,pos=4,labels="Top Non-Ivy+ CA Emp.",cex = .6,font=2,col="gray50")
    text(x=10,y=52,pos=4,labels="Ivy+ CA Emp.",cex = .6,font=2,col="blue")
    text(x=7.3,y=28,pos=4,labels="Ivy+ >$100k if Emp.",cex = .6,font=2,col="forestgreen")
    text(x=8,y=6,pos=4,labels="Top Non-Ivy+ >$100k if Emp.",cex = .6,font=2,col="red")
    legend(12.85,16,c("96-97 URM App.","98-99 URM App."),lty=c(1,2),cex=.6)
  #Same thing, but Black-only
  R<-data.frame(Time=6:16,Emp=NA,K100=NA)
  for(i in 6:16){
    R$Emp[i-5]<-mean(A[A$IvyEnrFirstIncCC==1&A$CETHNICA=="B"&A$YEARAPAY%in%1998:1999,paste0("wage_sum_",i,"_Employed")])*100
    R$EmpNon[i-5]<-mean(A[A$IvyEnrFirstIncCC==0&A$CETHNICA=="B"&A$YEARAPAY%in%1998:1999,paste0("wage_sum_",i,"_Employed")])*100
    R$EmpS4[i-5]<-mean(A[A$IvyEnrFirstIncCC==0&A$CETHNICA=="B"&A$YEARAPAY%in%1998:1999&NA_to_F(A$S_Q==4),paste0("wage_sum_",i,"_Employed")])*100
    R$Emp96[i-5]<-mean(A[A$IvyEnrFirstIncCC==1&A$CETHNICA=="B"&A$YEARAPAY%in%1996:1997,paste0("wage_sum_",i,"_Employed")])*100
    R$EmpNon96[i-5]<-mean(A[A$IvyEnrFirstIncCC==0&A$CETHNICA=="B"&A$YEARAPAY%in%1996:1997,paste0("wage_sum_",i,"_Employed")])*100
    R$EmpS496[i-5]<-mean(A[A$IvyEnrFirstIncCC==0&A$CETHNICA=="B"&A$YEARAPAY%in%1996:1997&NA_to_F(A$S_Q==4),paste0("wage_sum_",i,"_Employed")])*100
    R$K100[i-5]<-mean(A[A$IvyEnrFirstIncCC==1&A$CETHNICA=="B"&A$YEARAPAY%in%1998:1999,paste0("wage_sum_",i,"_100000")])*100
    R$K10096[i-5]<-mean(A[A$IvyEnrFirstIncCC==1&A$CETHNICA=="B"&A$YEARAPAY%in%1996:1997,paste0("wage_sum_",i,"_100000")])*100
    R$K100S4[i-5]<-mean(A[A$IvyEnrFirstIncCC==0&A$CETHNICA=="B"&A$YEARAPAY%in%1998:1999&NA_to_F(A$S_Q==4),paste0("wage_sum_",i,"_100000")])*100
    R$K100S496[i-5]<-mean(A[A$IvyEnrFirstIncCC==0&A$CETHNICA=="B"&A$YEARAPAY%in%1996:1997&NA_to_F(A$S_Q==4),paste0("wage_sum_",i,"_100000")])*100
  } 
  plot(R$Time,R$Emp,type="l",pch=NA,xlab="Years After Application",ylab="Proportion of Students",xlim=c(6,16),ylim=c(0,80),cex.lab=1.25,cex.axis=1.2,mgp=c(2.6,1,0),col="blue",lty=2) #Single-Diff
    lines(c(-10000,10000000),c(0,0))
    lines(R$Time,R$Emp96,col="blue")
    lines(R$Time,R$EmpNon,lty=2)
    lines(R$Time,R$EmpNon96)
    lines(R$Time,R$EmpS4,col="gray50",lty=2)
    lines(R$Time,R$EmpS496,col="gray50")
    lines(R$Time,R$K100/R$Emp*100,col="forestgreen",lwd=1.5,lty=2)
    lines(R$Time,R$K10096/R$Emp96*100,col="forestgreen",lwd=1.5)
    lines(R$Time,R$K100S4/R$EmpS4*100,col="red",lwd=1.5,lty=2)
    lines(R$Time,R$K100S496/R$EmpS496*100,col="red",lwd=1.5)
    text(x=12,y=69,pos=4,labels="Non-Ivy+ CA Emp.",cex = .6,font=2)
    text(x=5.5,y=58,pos=4,labels="Top Non-Ivy+ CA Emp.",cex = .45,font=2,col="gray50")
    text(x=10,y=50,pos=4,labels="Ivy+ CA Emp.",cex = .6,font=2,col="blue")
    text(x=7.3,y=28,pos=4,labels="Ivy+ >$100k if Emp.",cex = .6,font=2,col="forestgreen")
    text(x=8,y=6,pos=4,labels="Top Non-Ivy+ >$100k if Emp.",cex = .6,font=2,col="red")
    legend(12.85,16,c("96-97 Black App.","98-99 Black App."),lty=c(1,2),cex=.6)
  mean(A$IvyEnrFirstIncCC[A$URM&A$YEARAPAY%in%1998:1999],na.rm=T)-mean(A$IvyEnrFirstIncCC[A$URM&A$YEARAPAY%in%1996:1997],na.rm=T)
  mean(A$IvyEnrFirstIncCC[A$CETHNICA=="B"&A$YEARAPAY%in%1998:1999],na.rm=T)-mean(A$IvyEnrFirstIncCC[A$CETHNICA=="B"&A$YEARAPAY%in%1996:1997],na.rm=T)
  mean(R$K100/R$Emp*100) ; mean(R$K100S4/R$EmpS4*100)
  mean(R$Emp) ; mean(R$EmpS4)
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces the statistics presented in Table A-9.
  ##########################################################################################################
  ##########################################################################################################
  vs<-c("SAT","SATII_W","SATII_M","SATII_Other")
  for(v in vs){
    temp<-lm(as.formula(paste0(v,"~",paste0(vs[!vs%in%v],collapse="+"))),data=A[A$YEARAPAY%in%1996:1997,])
    A[,paste0(v,"_pred")]<-predict(temp,A)
    A[!is.na(A[,v]),paste0(v,"_pred")]<-A[!is.na(A[,v]),v]
    temp<-lm(as.formula(paste0(v,"~factor(CPREVSCH)+hsgpa")),data=A[A$YEARAPAY%in%1996:1997,])
    check<-A$CPREVSCH%in%gsub("^.*[)]","",names(temp$coefficients))
    A[is.na(A[,paste0(v,"_pred")])&check,paste0(v,"_pred")]<-predict(temp,A[is.na(A[,paste0(v,"_pred")])&check,])
  }
  A$S<-round((A$hsgpa_censor*1000+A$SATIV+A$SATIM+A$SATII_W+A$SATII_M+A$SATII_Other)/10)*10
  A$S_pred<-round((A$hsgpa_censor*1000+A$SAT_pred+A$SATII_W_pred+A$SATII_M_pred+A$SATII_Other_pred)/10)*10
  x1<-A[A$YEARAPAY%in%c(1995,1998:1999),]
  b<-unique(c(quantile(x1$S_pred[x1$URM], probs = seq(0, 1, by = 1/100),na.rm=T)))
  x1$Bin_pred<-as.integer(cut(x1$S_pred, breaks=b,labels=1:(length(b)-1), include.lowest=TRUE))
  hi_bin<-data.frame() #Recalculate, versus 94-95
  for(i in 1:100){
    x<-x1[x1$Bin_pred%in%(i-14):(i+14),]
    x$weight<- 1-(abs(x$Bin_pred-i)/15)
    y<-data.frame(Post=c(0,0,1,1),URM=c(0,1,0,1),Bin=rep(i,4))
    for(v in c(vars)) y[,v]<-c(weighted.mean(x[x$Post==0&x$URM==0,v],x$weight[x$Post==0&x$URM==0],na.rm=T),
                               weighted.mean(x[x$Post==0&x$URM==1,v],x$weight[x$Post==0&x$URM==1],na.rm=T),
                               weighted.mean(x[x$Post==1&x$URM==0,v],x$weight[x$Post==1&x$URM==0],na.rm=T),
                               weighted.mean(x[x$Post==1&x$URM==1,v],x$weight[x$Post==1&x$URM==1],na.rm=T))
    hi_bin<-rbind(hi_bin,y)
    print(i)
  }
  bs<-table(x1$Bin_pred[x1$URM&x1$YEARAPAY%in%1998:1999])/2
  load("Figures/AppDecline_AllUC_Ests.Rda") #Note: This is produced below
    Ests$App<-Ests$X2+Ests$X5 ; Ests<-Ests[,c("X1","App")]
  E<-data.frame(S_pred=seq(4000,8000,10),X1=(round((seq(4000,8000,10)-99)/200)*200))
  for(i in 1:401) E$Pop[i]<-mean(x1$S_pred[x1$URM&x1$YEARAPAY==1995]==E$S_pred[i],na.rm=T)
  E<-merge(E,x1[!duplicated(x1$S_pred),c("S_pred","Bin_pred")])
  camp<-data.frame()
  for(v1 in c("UCBEnrFirstIncCC","UCLAEnrFirstIncCC","UCSDEnrFirstIncCC","UCSBEnrFirstIncCC","UCIEnrFirstIncCC","UCDEnrFirstIncCC","UCSCEnrFirstIncCC","UCREnrFirstIncCC")){
    vs<-hi_bin[hi_bin$Post&hi_bin$URM,v1]-hi_bin[!hi_bin$Post&hi_bin$URM,v1]-(hi_bin[hi_bin$Post&!hi_bin$URM,v1]-hi_bin[!hi_bin$Post&!hi_bin$URM,v1])
    inc<-0 ; dec<-0
    for(i in 1:100) if(vs[i]<=0) dec<-dec-(vs[i]*bs[i]) else inc<-inc+(vs[i]*bs[i])
    app<-merge(E,data.frame(Bin_pred=1:100,
                            Change=0-hi_bin[!hi_bin$Post&hi_bin$URM,v1]-(hi_bin[hi_bin$Post&!hi_bin$URM,v1]-hi_bin[!hi_bin$Post&!hi_bin$URM,v1])))
    app<-ddply(app,.(X1),function(x) weighted.mean(x$Change,x$Pop))
    app<-merge(app,Ests)
    incapp<-0 ; decapp<-0
    for(i in 1:nrow(app)) if(app$App[i]<=0) decapp<-decapp+(app$V1[i]*app$App[i]) else incapp<-incapp-(app$V1[i]*app$App[i])
    camp<-rbind(camp,data.frame(C=v1,inc=inc,dec=dec,incapp=incapp,decapp=decapp))
    print(paste(v1,round(inc),round(dec),round(incapp),round(decapp)))
  }
  write.csv(camp,file="Figures/EnrChange_Total.csv")
  
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces Figure A-6.
  ##########################################################################################################
  ##########################################################################################################
  if(0==1){ #Only run once
    reg<-lm(logincome_parent~factor(CPREVSCH)+factor(PERM_ZIP)+SATIM_nom+SATIV_nom+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m+factor(SATII_1_nom)+factor(PAR_OCCUPN_KEY)+Female+factor(educ_lvl_max),A[A$YEARAPAY%in%1996:1997,])
    save(reg,file=paste0(secure_derived,"Income_Reg_AA.Rda"))
  } else load(paste0(secure_derived,"Income_Reg_AA.Rda"))
  cond<-rep(T,nrow(A)) ; for(var in c("CPREVSCH","PERM_ZIP","SATII_1_nom","PAR_OCCUPN_KEY")) cond<-cond&A[,var]%in%gsub("^.*[)]","",names(reg$coefficients)[grepl(var,names(reg$coefficients))])
  A$income_pred[cond]<-exp(predict(reg,A[cond,])) ; rm(reg) ; gc()  #Goes to 45GB
  A$Post<-A$YEARAPAY>1997
  A$income_topcode<-A$income_parent ; A$income_topcode[is.na(A$income_parent)]<-400000
  A$income_updated<-A$income_parent ; A$income_updated[is.na(A$income_parent)]<-A$income_pred[is.na(A$income_parent)]
  cor(A$income_parent,A$income_pred,use="complete.obs") #0.57
  
  hi<-data.frame()
  x<-A[A$YEARAPAY%in%1996:1999&(A$ENR01_Orig==1|A$ENR04_Orig==1),]
  for(i in seq(1000,300000,1000)){
    bw<-80000
    y<-data.frame(Post=c(0,0,1,1),URM=c(0,1,0,1),Bin=rep(i,4)/1000)
    y$income_parent<-c(mean(inrange(x$income_updated[x$Post==0],i-bw,i+bw-1)&x$URM[x$Post==0]==0,na.rm=T)*100,
                       mean(inrange(x$income_updated[x$Post==0],i-bw,i+bw-1)&x$URM[x$Post==0]==1,na.rm=T)*100,
                       mean(inrange(x$income_updated[x$Post==1],i-bw,i+bw-1)&x$URM[x$Post==1]==0,na.rm=T)*100,
                       mean(inrange(x$income_updated[x$Post==1],i-bw,i+bw-1)&x$URM[x$Post==1]==1,na.rm=T)*100)
    hi<-rbind(hi,y)
  }
  for(z in c(0,1)) hi$income_parent[hi$Post==z]<-hi$income_parent[hi$Post==z]/sum(hi$income_parent[hi$Post==z])*1000
  v<-"income_parent"
  tiff(paste0("Figures/IncDist.tiff"),5,4,"in",res=300)
    plot(hi$Bin[hi$Post&hi$URM],hi[!hi$Post&hi$URM,v],type="l",pch=NA,xlab="Family Income ($'000s)",xlim=c(0,300),ylab="% of Total Enr.",ylim=c(0,8),cex.lab=1.25,cex.axis=1.2,mgp=c(2.6,1,0)) #Single-Diff
    lines(c(-10000,10000000),c(0,0))
    lines(hi$Bin[hi$Post&hi$URM],hi[hi$Post&hi$URM,v],lty=2,col="black")
    lines(hi$Bin[hi$Post&hi$URM],hi[!hi$Post&!hi$URM,v],col="gray50")
    lines(hi$Bin[hi$Post&hi$URM],hi[hi$Post&!hi$URM,v],lty=2,col="gray50")
    legend(305,8,c("URM Pre-209","URM Post-209","Non-URM Pre-209","Non-URM Post-209"),lty=c(1,2,1,2),col=c("black","black","gray50","gray50"),xjust=1,yjust=1,cex=.85)
  dev.off()
  tiff(paste0("Figures/IncDistChange.tiff"),5,4,"in",res=300)
    plot(hi$Bin[hi$Post&hi$URM],hi[hi$Post&hi$URM,v]-hi[!hi$Post&hi$URM,v],type="l",pch=NA,xlab="Family Income ($'000s)",xlim=c(0,300),ylab="Change in % Total Enr.",ylim=c(-1.5,1),cex.lab=1.25,cex.axis=1.2,mgp=c(2.6,1,0)) #Single-Diff
    lines(c(-10000,10000000),c(0,0))
    lines(hi$Bin[hi$Post&!hi$URM],hi[hi$Post&!hi$URM,v]-hi[!hi$Post&!hi$URM,v],col="gray50")
    lines(hi$Bin[hi$Post&hi$URM],hi[hi$Post&!hi$URM,v]-hi[!hi$Post&!hi$URM,v]+(hi[hi$Post&hi$URM,v]-hi[!hi$Post&hi$URM,v]),lty=1,col="blue",lwd=2)
    abline(v=quantile(x$income_updated[x$YEARAPAY%in%1996:1997]/1000,c(.25),na.rm=T),col="gray50",lty=2)
    abline(v=quantile(x$income_updated[x$YEARAPAY%in%1996:1997]/1000,c(.50),na.rm=T),col="gray50",lty=2)
    abline(v=quantile(x$income_updated[x$YEARAPAY%in%1996:1997]/1000,c(.75),na.rm=T),col="gray50",lty=2)
    legend(305,-1.52,c("Change in URM","Change in Non-URM","Net Change"),col=c("black","gray50","blue"),xjust=1,yjust=0,lwd=c(1,1,2),cex=.85)
  dev.off()
  
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces the statistics presented in Tables I, A-3, A-4, and A-5.
  ##########################################################################################################
  ##########################################################################################################
  D<-data.frame(matrix("",32,80),stringsAsFactors = F)
  for(i in c(1,3:9)) A[,paste0("ENRR0",i)]<-A[,paste0("ENR0",i,"_Orig")]
  A$Pre<-NA ; A$Pre[A$YEARAPAY%in%1996:1997]<-1 ; A$Pre[A$YEARAPAY%in%1998:1999]<-0 ; A$Pre[A$YEARAPAY%in%1994:1995]<-2
  A$APP0100<-1 ; A$ADM0100<-A$UCAdm ; A$ENRR0100<-apply(A[,grepl("ENRR",names(A))],1,sum)>0
  
  r<-1 ; for(v in c("Num","SAT","GPA")){
    for(s in c(1,4,6,3,9,8,7,5,100)){
      c<-1 ; for(u in list(A$URM==1,A$URM==0,A$CETHNICA=="P",A$CETHNICA%in%c("D","E","G","H","L","N","V"),A$Black,A$Hispanic)){
        for(t in c("APP","ADM","ENRR")){
          for(p in c(2,1,0)){
            if(v=="Num"&t=="APP") D[r,c]<-ex(nrow(A[NA_to_F(u)&NA_to_F(A$Pre==p)&NA_to_F(A[,paste0(t,"0",s)]==1),])/2,0,comma=T)
            if(v=="Num"&t!="APP") D[r,c]<-ex(mean(A[NA_to_F(u)&NA_to_F(A$Pre==p)&NA_to_F(A[,paste0("APP0",s)]==1),paste0(t,"0",s)],na.rm=T)*100,1,comma=T)
            if(v=="SAT") D[r,c]<-ex(mean(A[NA_to_F(u)&NA_to_F(A$Pre==p)&NA_to_F(A[,paste0(t,"0",s)]==1),"SAT"],na.rm=T),0)
            if(v=="GPA") D[r,c]<-ex(mean(A[NA_to_F(u)&NA_to_F(A$Pre==p)&NA_to_F(A[,paste0(t,"0",s)]==1),"hsgpa"],na.rm=T),2)
            if(v=="HS") D[r,c]<-ex(length(unique(A[NA_to_F(u)&NA_to_F(A$Pre==p)&NA_to_F(A[,paste0(t,"0",s)]==1),"CPREVSCH"])),0)
            c<-c+1
          }
          c<-c+1
        }
        c<-c+1
      }
      r<-r+1
    }
    r<-r+2
  }
  write.csv(D,file="Figures/Prop209_Descriptives.csv")
  
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces the statistics presented in Table A-13.
  ##########################################################################################################
  ##########################################################################################################
  R<-data.frame(Major=gsub("NSC_Major_","",names(R5)),Change=unlist(R5[4,]),Change_sd=unlist(R5[5,]),val=as.numeric(gsub("[^0-9.-]","",unlist(R5[4,]))),stringsAsFactors = F)[-1,] ; R<-R[!grepl("Condit",R$Major),] ; R$Major<-gsub("([a-z])([A-Z])","\\1 \\2",R$Major) ; R<-R[order(R$val),]
  temp<-data.frame(Major=c(A$NSC_Major1[A$URM&A$YEARAPAY%in%1996:1997],A$NSC_Major[A$URM&A$YEARAPAY%in%1996:1997],A$NSC_Major3[A$URM&A$YEARAPAY%in%1996:1997]),Count=1) ; temp<-temp[!temp$Major%in%c("","NA"),] ; temp<-aggregate(Count~Major,temp,sum) ; majors<-read_csv("Data/Derived/NSC_Major_Codes_Categorized.csv") ; majors<-merge(majors[,!names(majors)=="Count"],temp) ; majors<-majors[order(majors$Count,decreasing = T),]
  R$Baseline<-"" ; R$Major1<-"" ; R$Major2<-"" ; R$Major3<-""
  for(v in row.names(R)){
    r<-which(row.names(R)==v)
    R$Baseline[r]<-ex(mean(A[A$YEARAPAY%in%1996:1997&A$URM,v],na.rm=T)*100,1)
    R[r,c("Major1","Major2","Major3")]<-titleCase(majors$Major[NA_to_F(majors$MajorSimp==R$Major[r])][1:3])
  }
  write.csv(R,file="Figures/Prop209_SpecificMajors.csv")
  
  
  #Test for "warming," as in Antonovics and Sander (2012)
  for(v in names(A)[grepl("^ADM0",names(A))]) A[,v]<-NA_to_F(A[,v])
  A$ADM_Set<-paste0(A$ADM01,A$ADM03,A$ADM04,A$ADM05,A$ADM06,A$ADM07,A$ADM07,A$ADM08,A$ADM09)
  summary(felm(UCEnr~URM*Prop209+SATIM+SATIV+hsgpa+SATIIW+SATII_M+SATII_M2_Indic+SATII_Other|CPREVSCH+factor(SATII_1)+ADM_Set,A[A$YEARAPAY%in%1996:1999,]),robust=T)
}











if(UC_Campus_Ests){
  load(paste0(secure_derived,"AffAct_Data.Rda"))
  A<-A[A$CQUARTERAP==2&A$SCHTYPE%in%c("A","B")&!is.na(A$hsgpa),]
  A$Prop209<-A$YEARAPAY>1997 ; A$Prop209[!inrange(A$YEARAPAY,1996,1999)]<-NA
  A$URM<-A$CETHNICA%in%c("A","B","C","J")
  A$APP99<-1
  A$ADM99<-NA_to_F(A$ADM01==1)|NA_to_F(A$ADM03==1)|NA_to_F(A$ADM04==1)|NA_to_F(A$ADM05==1)|NA_to_F(A$ADM06==1)|NA_to_F(A$ADM07==1)|NA_to_F(A$ADM08==1)|NA_to_F(A$ADM09==1)
  A$ENR99<-NA_to_F(A$ENR01==1)|NA_to_F(A$ENR03==1)|NA_to_F(A$ENR04==1)|NA_to_F(A$ENR05==1)|NA_to_F(A$ENR06==1)|NA_to_F(A$ENR07==1)|NA_to_F(A$ENR08==1)|NA_to_F(A$ENR09==1)
  A$URMPre209<-A$URM==1&A$Prop209==0&A$YEARAPAY!=1997
  A$URMPost209<-A$URM==1&A$Prop209==1
  A$Time<-0 ; A$Time[A$YEARAPAY==1997]<-1 ; A$Time[A$YEARAPAY>1997]<-2
  A$Asian<-A$CETHNICA%in%c("D","E","G","H","L","N","V")

  ##########################################################################################################
  ##########################################################################################################
  #The following code produces the statistics presented in Tables C-1, C-2, and A-6.
  ##########################################################################################################
  ##########################################################################################################
  
  R<-matrix("",nrow=55,ncol=9)
  c<-1
  for(i in c("01","04","06","08","09","03","07","05","99")){
    R<-genmat(R,1,c,paste0("APP",i),A,numdec=1)
    R<-genmat(R,11,c,paste0("ADM",i),A[A[,paste0("APP",i)]==1,],numdec=1)
    R<-genmat(R,21,c,paste0("ENR",i),A[A[,paste0("APP",i)]==1,],numdec=1)
    R<-genmat(R,31,c,paste0("ENR",i),A[NA_to_F(A[,paste0("ADM",i)]==1),],numdec=1)
    c<-c+1 ; print(i)
  }
  write.csv(R,file="Figures/Prop209_Effect.csv")
  
  #Now do the same comparing 95 to 98
  A$Prop209<-A$YEARAPAY>1997 ; A$Prop209[!inrange(A$YEARAPAY,1994,2001)]<-NA
  A$Prop209[!A$YEARAPAY%in%c(1994:1995,1998:1999)]<-NA
  R<-matrix("",nrow=55,ncol=9)
  c<-1
  for(i in c("01","04","06","08","09","03","07","05","99")){ #Ordered by 1998 overall admissions rate
    R<-genmat(R,1,c,paste0("APP",i),A,numdec=1)
    R<-genmat(R,11,c,paste0("ADM",i),A[A[,paste0("APP",i)]==1,],numdec=1)
    R<-genmat(R,21,c,paste0("ENR",i),A[A[,paste0("APP",i)]==1,],numdec=1)
    R<-genmat(R,31,c,paste0("ENR",i),A[NA_to_F(A[,paste0("ADM",i)]==1),],numdec=1)
    c<-c+1 ; print(i)
  }
  write.csv(R,file="Figures/Prop209_Effect_9495.csv")
  
  for(e in c("Black","Hispanic")){
    A$URM<-NA ; A$URM[!A$CETHNICA%in%c("A","B","C","J")]<-F
    if(e=="Black") A$URM[A$CETHNICA%in%c("B")]<-T
    if(e=="Hispanic") A$URM[A$CETHNICA%in%c("C","J")]<-T
    R<-matrix("",nrow=55,ncol=9)
    c<-1
    for(i in c("01","04","06","08","09","03","07","05","99")){
      R<-genmat(R,1,c,paste0("APP",i),A,numdec=1)
      R<-genmat(R,11,c,paste0("ADM",i),A[A[,paste0("APP",i)]==1,],numdec=1)
      R<-genmat(R,21,c,paste0("ENR",i),A[A[,paste0("APP",i)]==1,],numdec=1)
      R<-genmat(R,31,c,paste0("ENR",i),A[NA_to_F(A[,paste0("ADM",i)]==1),],numdec=1)
      c<-c+1 ; print(i)
    }
    write.csv(R,file=paste0("Figures/Prop209_Effect_",e,".csv"))
  }
}











if(STEM_Perform_Persist){
  load(paste0(secure_derived,"AffAct_Data.Rda"))
  load(paste0(secure_derived,"AffAct_Data_Courses.Rda"))

  A<-A[A$CQUARTERAP==2&A$SCHTYPE%in%c("A","B")&!is.na(A$hsgpa),]
  A$Prop209<-A$YEARAPAY>1997 ; A$Prop209[!inrange(A$YEARAPAY,1995,2000)]<-NA

  C$Term[grepl("Fall",C$Term,ignore.case=T)]<-"Fall"
    C$Term[grepl("Summer",C$Term,ignore.case=T)]<-"Summer"
    C$Term[grepl("Spring",C$Term,ignore.case=T)]<-"Spring"
    C$Term[grepl("Winter",C$Term,ignore.case=T)]<-"Winter"
  C$time<-C$Year-C$YEARAPAY+0.1*(C$Term=="Spring")+0.3*(C$Term=="Fall")
  C<-C[inrange(C$time,.3,2.3)&C$AL_KEY%in%A$AL_KEY,]
  
  #Identify science courses
  C$Chem1<-C$Course%in%c("1 CHEM 1A","3 Chemistry 002A","5 Chemistry 001A","8 CHEM       1A ","7 CHEM 1B")
  chem2<-c("1 CHEM 1B","3 Chemistry 002C","5 Chemistry 001B","8 CHEM       1B ","7 CHEM 1C")
    C$Chem2<-C$Course%in%c(chem2,"1 CHEM 3A","1 CHEM 12A") #Can just jump to Orgo at Berkeley
  chem3<-c("1 CHEM 3A","1 CHEM 112A","3 Chemistry 008A","3 Chemistry 118A","5 Chemistry 112A","8 CHEM     107A ","8 CHEM       6A ","7 CHEM 108A","7 CHEM 112A")
    C$Chem3<-C$Course%in%chem3
  chem4<-c("1 CHEM 3B","1 CHEM 112B","3 Chemistry 008B","3 Chemistry 118B","5 Chemistry 112B","8 CHEM     107B ","8 CHEM       6B ","7 CHEM 108B","7 CHEM 112B")
    C$Chem4<-C$Course%in%chem4
  C$Bio1<-C$Course%in%c("1 BIOLOGY 1B","3 Biological Sciences 001A","5 Biology 005A","8 MCDB       1A ","8 MCDB       4A ","8 BIOL       4A ","8 MCDB       5A ","7 BIOL 20A","7 BIOL 10") #Note: at Berkeley, 1B is taken before 1A by 73.3% of students who take both in this period
    C$Bio2<-C$Course%in%c("1 BIOLOGY 1A","3 Biological Sciences 001C","5 Biology 005C","8 EEMB       4C ","8 BIOL       4C ","8 EEMB       5C ","8 EEMB       2  ","7 BIOL 20C","7 BIOL 12")
  C$Math1<-C$Course%in%c("1 MATH 1A","1 MATH H1A","1 MATH 1AM","1 MATH 16A","3 Mathematics 016A","3 Mathematics 021A","5 Mathematics 009A","8 MATH       3A ","8 MATH      34A ","7 MATH 10A","7 MATH 11A","7 MATH 19A")
  C$Math2<-C$Course%in%c("1 MATH 1B","1 MATH H1B","1 MATH 1BM","1 MATH 16B","3 Mathematics 016B","3 Mathematics 021B","5 Mathematics 009B","8 MATH       3B ","8 MATH      34B ","7 MATH 10B","7 MATH 11B","7 MATH 19B")
  C$Math3<-C$Course%in%c("1 MATH 53","1 MATH H53","1 MATH 16B","3 Mathematics 016C","3 Mathematics 021C","5 Mathematics 009C","8 MATH       3C ","7 MATH 11C","7 MATH 22","7 MATH 23A")
  C$Physics1<-C$Course%in%c("1 PHYSICS 8A","3 Physics 001A","3 Physics 005A","3 Physics 007A","3 Physics 009A","3 Physics 009HA","5 Physics 002A","8 PHYS       6A ","7 PHYS 6A","7 PHYS 5A","7 PHYS 7A")
    C$Physics2<-C$Course%in%c("1 PHYSICS 8B","3 Physics 001B","3 Physics 005C","3 Physics 007C","3 Physics 009C","3 Physics 009HC","5 Physics 002C","8 PHYS       6C ","7 PHYS 6C","7 PHYS 5C","7 PHYS 7B")
  compsci1<-c("1 COMPSCI 61A","3 Engineering Computer Science 020","3 Engineering Computer Science 030","5 Computer Science 010","8 CMPSC     10  ","7 CMPS 12A")
    C$CompSci1<-C$Course%in%compsci1
  compsci2<-c("1 COMPSCI 61B","3 Engineering Computer Science 040","5 Computer Science 012","8 CMPSC     20  ","7 CMPS 12B")
    C$CompSci2<-C$Course%in%compsci2
  compsci3<-c("1 COMPSCI 61C","3 Engineering Computer Science 050","3 Engineering Electrical & Compu 070","5 Computer Science 014","8 CMPSC     30  ","7 CMPE 12C","7 CMPS 101")
    C$CompSci3<-C$Course%in%compsci3

  #Prepare A
  A$URMPre209<-A$URM==1&A$Prop209==0 
  A$URMPost209<-A$URM==1&A$Prop209==1
  A$Black<-A$CETHNICA%in%c("B") ; A$Black[A$CETHNICA%in%c("A")]<-NA
  A$Hispanic<-A$CETHNICA%in%c("C","J") ; A$Hispanic[A$CETHNICA%in%c("A")]<-NA
  
  #Define a UC-only STEM variable
  stem<-unlist(read.csv("Data/Raw/Definition of STEM.csv",stringsAsFactors = F)) ; stem<-gsub("^.*[^0-9]([0-9]{6})[^0-9].*$","\\1",gsub("[.]","",stem[grepl("[.]",stem)])) ; stem[nchar(stem)==5]<-paste0("0",stem[nchar(stem)==5]) ;  for(z in stem[grepl("X",stem)]) stem<-c(stem,paste0(gsub("X","",z),addzeros(1:9999,4)))
  A$STEM_NoNSC<-NA_to_F(A$CIP6RECENT%in%stem)*100 ; A$STEM_NoNSC[is.na(A$HASGRAD6)]<-NA
  
  C<-C[order(C$AL_KEY,C$time),]
  
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces Figures VIII and A-17 and the statistics presented in Tables G-1, G-2, A-16, and A-1.
  ##########################################################################################################
  ##########################################################################################################
  Results<-data.frame(matrix("",nrow=80,ncol=36),stringsAsFactors = F)
  col<-1 ; for(version in c("TakeSTEM","TakeSTEM_All")){
    Aggf<-data.frame() ; Aggi<-data.frame() ; Aggg<-data.frame()
    data<-A[A$AL_KEY%in%C$AL_KEY,]
    if(version=="TakeSTEM") data<-data[data$ENR01_Orig==1,]
    if(version=="TakeSTEM_All") data<-data
    c<-merge(C,data[,c("AL_KEY","Prop209","URMPost209","hsgpa","SATIM_m","SATIV_m","SATIIW_m","SATII_M_m","SATII_Other_m","SATIM_nom","SATIV_nom","SATIIW_nom","SATII_M_nom","SATII_M2_Indic_nom","SATII_Other_nom","SATII_1_nom","Female","HSCRSE_TYP_ACADHONRS_12","CPREVSCH","Major01","CMAJPRP01","STEM_NSCOnly","STEM_NoNSC","CETHNICA","Black","Hispanic")]) ; c$Course_FE<-paste0(c$Course,c$Year,c$Term) ; c<-c[order(c$AL_KEY,c$time),]
    for(v in c("Chem1","Chem2","Chem3","Chem4","Bio1","Bio2","Physics1","Physics2","CompSci1","CompSci2","CompSci3")){
      names(Results)[col]<-v
      data$Indic<-(data$AL_KEY%in%(c$AL_KEY[(c[,v]==1)&c$Pass]))*100
      vbefore<-v ; for(i in 2:4) vbefore<-gsub(i,i-1,vbefore)
      if(!grepl(1,v)){
        temp<-c[c[,vbefore]==1,] ; temp<-temp[!duplicated(temp$AL_KEY),c("AL_KEY","Course","Course_FE")]
        data<-data[,!names(data)%in%c("Course","Course_FE")] ; data<-merge(data,temp,all.x=T)
      }
      if(!grepl(1,v)) cond<-data$AL_KEY%in%c$AL_KEY[c[,vbefore]] else cond<-rep(T,nrow(data))
      ycond<-data$YEARAPAY%in%1996:1999
      t1<-summary(felm(Indic~URM+URMPost209 | YEARAPAY,data=data[cond&ycond,],exactDOF=T),robust=T) ; t<-t1$coefficients ; t<-t[c("URMTRUE","URMPost209TRUE"),]
      r<-1 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(data$Indic[cond&ycond],na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T)
      t1<-summary(felm(Indic~URM+URMPost209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY,data=data[cond&ycond,],exactDOF=T),robust=T) ; t<-t1$coefficients ; t<-t[c("URMTRUE","URMPost209TRUE"),]
      r<-10 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(data$Indic[cond&ycond],na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T)
      
      c_data<-c[c[,v]==1,] ; c_data<-c_data[!duplicated(c_data$AL_KEY),]
      cycond<-c_data$YEARAPAY%in%1996:1999
      t1<-summary(felm(GPA~URM+URMPost209 | YEARAPAY,data=c_data[cycond,],exactDOF=T),robust=T) ; t<-t1$coefficients ; t<-t[c("URMTRUE","URMPost209TRUE"),]
      r<-20 ; Results[r,col]<-ex(t[1,1],2) ; Results[r+1,col]<-ex(t[1,2],2,par=T) ; Results[r+3,col]<-ex(t[2,1],2) ; Results[r+4,col]<-ex(t[2,2],2,par=T) ; Results[r+7,col]<-ex(mean(c_data[cycond,]$GPA,na.rm=T),2) ; Results[r+6,col]<-ex(t1$N,0,comma=T)
      t1<-summary(felm(GPA~URM+Prop209+URMPost209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY,data=c_data[cycond,],exactDOF=T),robust=T) ; t<-t1$coefficients ; t<-t[c("URMTRUE","URMPost209TRUE"),]
      r<-30 ; Results[r,col]<-ex(t[1,1],2) ; Results[r+1,col]<-ex(t[1,2],2,par=T) ; Results[r+3,col]<-ex(t[2,1],2) ; Results[r+4,col]<-ex(t[2,2],2,par=T) ; Results[r+7,col]<-ex(mean(c_data[cycond,]$GPA,na.rm=T),2) ; Results[r+6,col]<-ex(t1$N,0,comma=T)
      
      col<-col+1
      
      tempi<-data[cond,] ; tempi$Class<-v
      if(grepl("1",v)){
        tempi<-data.frame()
      }
      tempg<-c_data ; tempg$Class<-v
      Aggi<-rbind(Aggi,tempi) ; Aggg<-rbind(Aggg,tempg)
    }
    Aggf<-data
    Aggf$Indic<-(data$AL_KEY%in%(c$AL_KEY[(c$Chem1|c$Chem2|c$Chem3|c$Chem4|c$Bio1|c$Bio2|c$Physics1|c$Physics2|c$CompSci1|c$CompSci2|c$CompSci3)&!is.na(c$GPA)]))*100
    Aggg$SAT<-Aggg$SATIM_nom+Aggg$SATIV_nom ; Aggg$SAT[Aggg$SAT==0]<-NA ; Aggg$Count<-1
    Aggg<-ddply(Aggg,.(Course_FE),function(x){
      for(i in 1:nrow(x)) x$SAT_Pctile[i]<-mean(x$SAT[i]>x$SAT[-i],na.rm=T)*100
      return(x) 
    })
    #Now do the aggregate
    for(d in c("Aggi","Aggg","Aggf")){
      data<-get(d)
      data$Count<-1
      temp<-aggregate(Count~AL_KEY,data,sum) ; temp$Count<-1/temp$Count ; names(temp)[2]<-"Weight" 
      data<-merge(data[,!names(data)=="Weight"],temp)
      for(y in 1994:2002){
        data[,paste0("URM_",y)]<-data$URM*data$YEARAPAY==y
        data[,paste0("Black_",y)]<-data$Black*data$YEARAPAY==y
        data[,paste0("Hispanic_",y)]<-data$Hispanic*data$YEARAPAY==y
      } 
      data$BlackPost209<-data$Black&data$Prop209
      data$HispanicPost209<-data$Hispanic&data$Prop209
      assign(paste0(d,"_orig"),data) ; data<-data[data$YEARAPAY%in%1996:1999,]
      assign(d,data)
    }
    names(Results)[col]<-"ALL"
    
    t1<-summary(felm(Indic~URM+URMPost209 | YEARAPAY|0|AL_KEY+Course,data=Aggi,exactDOF=T),robust=T,weight=Aggi$Weight) ; t<-t1$coefficients ; t<-t[c("URMTRUE","URMPost209TRUE"),]
    r<-1 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(Aggi$Indic,na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T)
    t1<-summary(felm(Indic~URM+URMPost209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY|0|AL_KEY+Course,data=Aggi,exactDOF=T),robust=T,weight=Aggi$Weight) ; t<-t1$coefficients ; t<-t[c("URMTRUE","URMPost209TRUE"),]
    r<-10 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(Aggi$Indic,na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T)
    
    t1<-summary(felm(GPA~URM+URMPost209 | YEARAPAY|0|AL_KEY+Course,data=Aggg,exactDOF=T),robust=T,weight=Aggg$Weight) ; t<-t1$coefficients ; t<-t[c("URMTRUE","URMPost209TRUE"),]
    r<-20 ; Results[r,col]<-ex(t[1,1],2) ; Results[r+1,col]<-ex(t[1,2],2,par=T) ; Results[r+3,col]<-ex(t[2,1],2) ; Results[r+4,col]<-ex(t[2,2],2,par=T) ; Results[r+7,col]<-ex(mean(Aggg$GPA,na.rm=T),2) ; Results[r+6,col]<-ex(t1$N,0,comma=T)
    t1<-summary(felm(GPA~URM+Prop209+URMPost209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY|0|AL_KEY+Course,data=Aggg,exactDOF=T,weight=Aggg$Weight),robust=T) ; t<-t1$coefficients ; t<-t[c("URMTRUE","URMPost209TRUE"),] 
    r<-30 ; Results[r,col]<-ex(t[1,1],2) ; Results[r+1,col]<-ex(t[1,2],2,par=T) ; Results[r+3,col]<-ex(t[2,1],2) ; Results[r+4,col]<-ex(t[2,2],2,par=T) ; Results[r+7,col]<-ex(mean(Aggg$GPA,na.rm=T),2) ; Results[r+6,col]<-ex(t1$N,0,comma=T)
    
    
    t1<-summary(felm(SAT_Pctile~URM+URMPost209 | YEARAPAY|0|AL_KEY+Course,data=Aggg,exactDOF=T),robust=T,weight=Aggg$Weight) ; t<-t1$coefficients ; t<-t[c("URMTRUE","URMPost209TRUE"),]
    r<-40 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(Aggg$SAT_Pctile,na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T)
    t1<-summary(felm(SAT_Pctile~URM+Prop209+URMPost209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY|0|AL_KEY+Course,data=Aggg,exactDOF=T,weight=Aggg$Weight),robust=T) ; t<-t1$coefficients ; t<-t[c("URMTRUE","URMPost209TRUE"),]
    r<-50 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(Aggg$SAT_Pctile,na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T)
    
    t1<-summary(felm(NA_to_F(STEM_NoNSC)~URM+URMPost209 | YEARAPAY,data=Aggf,exactDOF=T),robust=T,weight=Aggf$Weight) ; t<-t1$coefficients ; t<-t[c("URMTRUE","URMPost209TRUE"),]
    r<-60 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(Aggf$STEM_NoNSC,na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T)
    t1<-summary(felm(NA_to_F(STEM_NoNSC)~URM+URMPost209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY,data=Aggf,exactDOF=T),robust=T,weight=Aggf$Weight) ; t<-t1$coefficients ; t<-t[c("URMTRUE","URMPost209TRUE"),]
    r<-70 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(Aggf$STEM_NoNSC,na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T)
    
    if(version=="TakeSTEM_All"){
      #Now add results for appendix table
      col<-col+1
      check<-Aggi$YEARAPAY%in%1996:1999
      gpa_med<-quantile(Aggi$hsgpa[Aggi$Campus_ID==1&Aggi$URM&check],2/3)
      sat_med<-quantile((Aggi$AIS-Aggi$hsgpa_censor*1000)[Aggi$Campus_ID==1&Aggi$URM&check],2/3,na.rm=T)
      t1<-summary(felm(Indic~URM+URMPost209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY|0|AL_KEY+Course,data=Aggi[Aggi$Campus_ID==7&check,],exactDOF=T,weight=Aggi$Weight[Aggi$Campus_ID==7&check]),robust=T) ; t<-t1$coefficients[c("URMTRUE","URMPost209TRUE"),] #UC Santa Barbara
      r<-1 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(Aggi$Indic[Aggi$Campus_ID==7&check],na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T) ; col<-col+1
      t1<-summary(felm(Indic~URM+URMPost209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY|0|AL_KEY+Course,data=Aggi[Aggi$Campus_ID==3&check,],exactDOF=T,weight=Aggi$Weight[Aggi$Campus_ID==3&check]),robust=T) ; t<-t1$coefficients[c("URMTRUE","URMPost209TRUE"),] #UC Davis
      r<-1 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(Aggi$Indic[Aggi$Campus_ID==3&check],na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T) ; col<-col+1
      t1<-summary(felm(Indic~URM+URMPost209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY|0|AL_KEY+Course,data=Aggi[Aggi$Campus_ID==8&check,],exactDOF=T,weight=Aggi$Weight[Aggi$Campus_ID==8&check]),robust=T) ; t<-t1$coefficients[c("URMTRUE","URMPost209TRUE"),] #UC Santa Cruz
      r<-1 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(Aggi$Indic[Aggi$Campus_ID==8&check],na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T) ; col<-col+1
      t1<-summary(felm(Indic~URM+URMPost209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY|0|AL_KEY+Course,data=Aggi[Aggi$Campus_ID==5&check,],exactDOF=T,weight=Aggi$Weight[Aggi$Campus_ID==5&check]),robust=T) ; t<-t1$coefficients[c("URMTRUE","URMPost209TRUE"),] #UC Riverside
      r<-1 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(Aggi$Indic[Aggi$Campus_ID==5&check],na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T) ; col<-col+1
      t1<-summary(felm(Indic~URM+URMPost209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m+income_parent_zeros+income_Missing | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY+PAR_OCCUPN_KEY+educ_lvl_max|0,data=Aggi[Aggi$Campus_ID==1&check,],exactDOF=T,weight=Aggi$Weight[Aggi$Campus_ID==1&check]),robust=T) ; t<-t1$coefficients[c("URMTRUE","URMPost209TRUE"),] #Extra controls
      r<-1 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(Aggi$Indic[Aggi$Campus_ID==1&check],na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T) ; col<-col+1
      cond<-check&NA_to_F(Aggi$Campus_ID==1&(Aggi$AIS-Aggi$hsgpa_censor*1000)>=sat_med&Aggi$hsgpa>=gpa_med)
      t1<-summary(felm(Indic~URM+URMPost209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY|0|AL_KEY+Course,data=Aggi[cond,],exactDOF=T,weight=Aggi$Weight[cond]),robust=T) ; t<-t1$coefficients[c("URMTRUE","URMPost209TRUE"),]
      r<-1 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(Aggi$Indic[cond],na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T) ; col<-col+1
      cond<-check&NA_to_F(Aggi$Campus_ID==1&(Aggi$AIS-Aggi$hsgpa_censor*1000)>=sat_med&Aggi$hsgpa<gpa_med)
      t1<-summary(felm(Indic~URM+URMPost209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY|0|AL_KEY+Course,data=Aggi[cond,],exactDOF=T,weight=Aggi$Weight[cond]),robust=T) ; t<-t1$coefficients[c("URMTRUE","URMPost209TRUE"),]
      r<-1 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(Aggi$Indic[cond],na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T) ; col<-col+1
      cond<-check&NA_to_F(Aggi$Campus_ID==1&Aggi$hsgpa>=gpa_med&(Aggi$AIS-Aggi$hsgpa_censor*1000)<sat_med)
      t1<-summary(felm(Indic~URM+URMPost209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY|0|AL_KEY+Course,data=Aggi[cond,],exactDOF=T,weight=Aggi$Weight[cond]),robust=T) ; t<-t1$coefficients[c("URMTRUE","URMPost209TRUE"),]
      r<-1 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(Aggi$Indic[cond],na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T) ; col<-col+1
      cond<-check&NA_to_F(Aggi$Campus_ID==1&Aggi$hsgpa<gpa_med&(Aggi$AIS-Aggi$hsgpa_censor*1000)<sat_med)
      t1<-summary(felm(Indic~URM+URMPost209 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY|0|AL_KEY+Course,data=Aggi[cond,],exactDOF=T,weight=Aggi$Weight[cond]),robust=T) ; t<-t1$coefficients[c("URMTRUE","URMPost209TRUE"),]
      r<-1 ; Results[r,col]<-ex(t[1,1],1) ; Results[r+1,col]<-ex(t[1,2],1,par=T) ; Results[r+3,col]<-ex(t[2,1],1) ; Results[r+4,col]<-ex(t[2,2],1,par=T) ; Results[r+7,col]<-ex(mean(Aggi$Indic[cond],na.rm=T),1) ; Results[r+6,col]<-ex(t1$N,0,comma=T) ; col<-col+1
      
      #Now make some plots to show changes over time
      Aggi<-Aggi_orig ; Aggg<-Aggg_orig ; Aggf<-Aggf_orig
      
      select_all<-rbind(data.frame(summary(felm(SAT_Pctile~URM_1994+URM_1995+URM_1996+URM_1998+URM_1999+URM_2000+URM_2001+URM_2002 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY+CETHNICA|0|AL_KEY+Course,data=Aggg,exactDOF=T,weight=Aggg$Weight),robust=T)$coefficients[1:8,1:2]),data.frame(Estimate=0,Cluster.s.e.=0)) ; select_all<-select_all[order(c(1994:1996,1998:2002,1997)),]
      select_all_urm<-rbind(data.frame(summary(felm(SAT_Pctile~Black_1994+Black_1995+Black_1996+Black_1998+Black_1999+Black_2000+Black_2001+Black_2002+Hispanic_1994+Hispanic_1995+Hispanic_1996+Hispanic_1998+Hispanic_1999+Hispanic_2000+Hispanic_2001+Hispanic_2002 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY+CETHNICA|0|AL_KEY+Course,data=Aggg,exactDOF=T,weight=Aggg$Weight),robust=T)$coefficients[1:16,1:2]),data.frame(Estimate=c(0,0),Cluster.s.e.=c(0,0))) ; select_all_urm$Eth<-c(rep("Black",8),rep("Hispanic",8),"Black","Hispanic") ; select_all_urm<-select_all_urm[order(select_all_urm$Eth,c(1994:1996,1998:2002,1994:1996,1998:2002,1997,1997)),]
      
      indic_all<-rbind(data.frame(summary(felm(Indic~URM_1994+URM_1995+URM_1996+URM_1998+URM_1999+URM_2000+URM_2001+URM_2002 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY+CETHNICA|0|AL_KEY+Course,data=Aggi,exactDOF=T,weight=Aggi$Weight),robust=T)$coefficients[1:8,1:2]),data.frame(Estimate=0,Cluster.s.e.=0)) ; indic_all<-indic_all[order(c(1994:1996,1998:2002,1997)),]
      indic_all_urm<-rbind(data.frame(summary(felm(Indic~Black_1994+Black_1995+Black_1996+Black_1998+Black_1999+Black_2000+Black_2001+Black_2002+Hispanic_1994+Hispanic_1995+Hispanic_1996+Hispanic_1998+Hispanic_1999+Hispanic_2000+Hispanic_2001+Hispanic_2002 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY+CETHNICA|0|AL_KEY+Course,data=Aggi,exactDOF=T,weight=Aggi$Weight),robust=T)$coefficients[1:16,1:2]),data.frame(Estimate=c(0,0),Cluster.s.e.=c(0,0))) ; indic_all_urm$Eth<-c(rep("Black",8),rep("Hispanic",8),"Black","Hispanic") ; indic_all_urm<-indic_all_urm[order(indic_all_urm$Eth,c(1994:1996,1998:2002,1994:1996,1998:2002,1997,1997)),]
      
      grade_all<-rbind(data.frame(summary(felm(GPA~URM_1994+URM_1995+URM_1996+URM_1998+URM_1999+URM_2000+URM_2001+URM_2002 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY+CETHNICA|0|AL_KEY+Course,data=Aggg,exactDOF=T,weight=Aggg$Weight),robust=T)$coefficients[1:8,1:2]),data.frame(Estimate=0,Cluster.s.e.=0)) ; grade_all<-grade_all[order(c(1994:1996,1998:2002,1997)),]
      grade_all_urm<-rbind(data.frame(summary(felm(GPA~Black_1994+Black_1995+Black_1996+Black_1998+Black_1999+Black_2000+Black_2001+Black_2002+Hispanic_1994+Hispanic_1995+Hispanic_1996+Hispanic_1998+Hispanic_1999+Hispanic_2000+Hispanic_2001+Hispanic_2002 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY+CETHNICA|0|AL_KEY+Course,data=Aggg,exactDOF=T,weight=Aggg$Weight),robust=T)$coefficients[1:16,1:2]),data.frame(Estimate=c(0,0),Cluster.s.e.=c(0,0))) ; grade_all_urm$Eth<-c(rep("Black",8),rep("Hispanic",8),"Black","Hispanic") ; grade_all_urm<-grade_all_urm[order(grade_all_urm$Eth,c(1994:1996,1998:2002,1994:1996,1998:2002,1997,1997)),]
      deg_all<-rbind(data.frame(summary(felm(STEM_NoNSC~URM_1994+URM_1995+URM_1996+URM_1998+URM_1999+URM_2000+URM_2001+URM_2002 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY+CETHNICA,data=Aggf,exactDOF=T,weight=Aggf$Weight),robust=T)$coefficients[1:8,1:2]),data.frame(Estimate=0,Robust.s.e=0)) ; deg_all<-deg_all[order(c(1994:1996,1998:2002,1997)),]
      deg_all_urm<-rbind(data.frame(summary(felm(STEM_NoNSC~Black_1994+Black_1995+Black_1996+Black_1998+Black_1999+Black_2000+Black_2001+Black_2002+Hispanic_1994+Hispanic_1995+Hispanic_1996+Hispanic_1998+Hispanic_1999+Hispanic_2000+Hispanic_2001+Hispanic_2002 + SATIM_nom+SATIV_nom+hsgpa+SATIIW_nom+SATII_M_nom+factor(SATII_M2_Indic_nom)+SATII_Other_nom+SATIM_m+SATIV_m+SATIIW_m+SATII_M_m+SATII_Other_m | factor(CPREVSCH)+factor(SATII_1_nom)+YEARAPAY+CETHNICA,data=Aggf,exactDOF=T,weight=Aggf$Weight),robust=T)$coefficients[1:16,1:2]),data.frame(Estimate=c(0,0),Robust.s.e=c(0,0))) ; deg_all_urm$Eth<-c(rep("Black",8),rep("Hispanic",8),"Black","Hispanic") ; deg_all_urm<-deg_all_urm[order(deg_all_urm$Eth,c(1994:1996,1998:2002,1994:1996,1998:2002,1997,1997)),]
      
      for(v in list(
        list("select_all","Percentiles",c(-3,10),seq(-2.5,10,2.5),c("",0,"",5,"",10),c(-3,13),seq(-2.5,12.5,2.5),c(NA,0,NA,5,NA,10,NA)),
        list("indic_all","Percentage Points",c(-14,14),seq(-10,10,5),c(-10,"",0,"",10),c(-17,14),seq(-15,10,5),c(NA,-10,"",0,"",10)),
        list("grade_all","Grade Points",c(-0.18,0.15),seq(-.15,.15,.05),c(NA,"-0.1",NA,0,NA,"0.1",NA),c(-0.31,0.25),seq(-.3,.2,.1),c("-0.3",NA,"-0.1",NA,"0.1",NA)),
        list("deg_all","Percentage Points",c(-10,10),seq(-10,10,5),c(-10,-5,0,5,10),c(-10,10),seq(-10,10,5),c(-10,-5,0,5,10)))){
        tiff(paste0("Figures/AA_CourseEffects_",v[[1]],".tiff"),5,3.5,"in",res=300)
        plot(1994:2002, get(v[[1]])[,1],
             ylim=v[[3]],
             pch=16, xlab="", ylab=v[[2]],cex.lab=1.1,cex.axis=1.2,mgp=c(2.4,1,0),yaxt="n"
        )
        axis(2, labels =v[[5]] , at = v[[4]],cex.axis=1.05)
        lines(c(1997.5,1997.5),c(-100000,100000),col="gray60")
        lines(c(0,2100),rep(0,2),col="gray")
        points(1994:2002, get(v[[1]])[,1],pch=16)
        arrows(1994:2002,  get(v[[1]])[,1]- get(v[[1]])[,2]*1.96, 1994:2002,  get(v[[1]])[,1]+get(v[[1]])[,2]*1.96, length=0.05, angle=90, code=3)
        dev.off()
        o<-get(paste0(v[[1]],"_urm"))
        tiff(paste0("Figures/AA_CourseEffects_",v[[1]],"_MultiURM.tiff"),5,3.5,"in",res=300)
        plot(1994:2002-0.1, o[o$Eth=="Black",1],
             ylim=v[[6]],xlim=c(1993.9,2002.1),
             pch=17, xlab="", ylab=v[[2]],cex.lab=1.1,cex.axis=1.2,mgp=c(2.4,1,0),yaxt="n"
        )
        axis(2, labels =v[[8]] , at = v[[7]],cex.axis=1.05)
        lines(c(1997.5,1997.5),c(-100000,100000),col="gray60")
        lines(c(0,2100),rep(0,2),col="gray")
        points(1994:2002-.1, o[o$Eth=="Black",1],pch=17)
        arrows(1994:2002-.1, o[o$Eth=="Black",1]-o[o$Eth=="Black",2]*1.96, 1994:2002-.1, o[o$Eth=="Black",1]+o[o$Eth=="Black",2]*1.96, length=0.05, angle=90, code=3)
        points(1994:2002+.1, o[o$Eth=="Hispanic",1],pch=15,col="gray50")
        arrows(1994:2002+.1, o[o$Eth=="Hispanic",1]-o[o$Eth=="Hispanic",2]*1.96, 1994:2002+.1, o[o$Eth=="Hispanic",1]+o[o$Eth=="Hispanic",2]*1.96, length=0.05, angle=90, code=3,col="gray50")
        dev.off()
      }
    }
    col<-col+2
  }
  write.csv(Results,"Figures/Prop209_Courses_Sci_Berkeley.csv")
}
























if(App_Charts_and_Discontinuities){
  load(paste0(secure_derived,"AffAct_Data.Rda"))
    A$UnivGrad<-(A$UnivGrad_NSCOnly==1)|NA_to_F(A$HASGRAD5==1)
    A$UnivGrad_SixYears<-(A$UnivGrad_SixYears_NSCOnly==1)|NA_to_F(A$HASGRAD6==1)  
    
  #Identify 'normal' UC applicants to each campus; excludes colleges that ran separate admissions (following archival documentation)
  maj<-data.frame(read_excel("Data/Raw/Majors_to_CIP.xlsx",col_names = F),stringsAsFactors = F)[,c(3,5,6,7)]
    names(maj)<-c("Campus","School","Major","Major_Name")
    for(zzz in 1:2) maj$Major[nchar(maj$Major)<3]<-paste0(0,maj$Major[nchar(maj$Major)<3])
  for(s in c("01","03","04","05","06","07","08","09")) for(zzz in 1:2) A[NA_to_F(nchar(A[,paste0("CMAJPRP",s)])<3),paste0("CMAJPRP",s)]<-paste0(0,A[NA_to_F(nchar(A[,paste0("CMAJPRP",s)])<3),paste0("CMAJPRP",s)])
  A$LS01<-!A$CMAJPRP01%in%maj$Major[maj$Campus=="BK"&!maj$School%in%c("COLLEGE OF LETTERS & SCIENCE","GRADUATE DIVISION","UNAFFILIATED")]
  A$LS03<-!A$CMAJPRP03%in%maj$Major[maj$Campus=="DV"&!maj$School%in%c("GRADUATE DIVISION","COLLEGE OF LETTERS & SCIENCE","COLLEGE OF AGRIC & ENVIRON SCIENCES")]
  A$LS09<-!A$CMAJPRP09%in%maj$Major[maj$Campus=="IR"&!maj$School%in%c("GRADUATE DIVISION","SCHOOL OF HUMANITIES","SCHOOL OF THE ARTS","SCHOOL OF SOCIAL SCIENCES","SCHOOL OF PHYSICAL SCIENCES")]
  A$LS04<-!A$CMAJPRP04%in%maj$Major[maj$Campus=="LA"&!maj$School%in%c("COLLEGE OF LETTERS & SCIENCE","GRADUATE DIVISION")]
  A$LS05<-!A$CMAJPRP05%in%maj$Major[maj$Campus=="RV"&!maj$School%in%c("COLLEGE OF HUMANITIES & SOCIAL SCI","COLLEGE OF NATURAL & AGRICULTURAL SCI","GRADUATE DIVISION")]
  A$LS06<-T
  A$LS08<-!A$CMAJPRP08%in%maj$Major[maj$Campus=="SB"&!maj$School%in%c("COLLEGE OF LETTERS & SCIENCE","GRADUATE DIVISION")]
  A$LS07<-T
  
  for(v in names(A)[grepl("^wage_sum_",names(A))]){
    A[,v]<-winsor(A[,v],0.01) #Winsorize 1%
    A[,paste0(v,"_Uncondit")]<-A[,v] ; A[is.na(A[,paste0(v,"_Uncondit")]),paste0(v,"_Uncondit")]<-0
    A[,paste0(v,"_Employed")]<-NA_to_F(A[,v]>0)
    A[,paste0(v,"_75000")]<-NA_to_F(A[,v]>=75000)
    A[,paste0(v,"_100000")]<-NA_to_F(A[,v]>=100000)
    A[,paste0(v,"_150000")]<-NA_to_F(A[,v]>=150000)
    A[NA_to_F(A[,v]!=0),paste0(v,"_log")]<-log(A[NA_to_F(A[,v]!=0),v])
    A[,paste0(v,"_logUncond")]<-A[,paste0(v,"_log")] ; A[is.na(A[,paste0(v,"_logUncond")]),paste0(v,"_logUncond")]<-0
  }
  A$Total_Wages_NumYears<-apply(A[,grepl("wage_sum_([6-9]|1[0-6])_Uncondit",names(A))],1,function(x) sum(x>0))
  A$Total_Wages_logUncond<-log(1+apply(A[,grepl("wage_sum_([6-9]|1[0-6])_Uncondit",names(A))],1,function(x) sum(x))/14)
  A$Total_Wages_logCondit<-log(1+apply(A[,grepl("wage_sum_([6-9]|1[0-6])_Uncondit",names(A))],1,function(x) sum(x))/A$Total_Wages_NumYears)
  A$Total_Wages_Num75000<-apply(A[,grepl("wage_sum_([6-9]|1[0-6])_Uncondit",names(A))],1,function(x) sum(x>75000))
  A$Total_Wages_Num150000<-apply(A[,grepl("wage_sum_([6-9]|1[0-6])_Uncondit",names(A))],1,function(x) sum(x>150000))
  A$Total_Wages30s_NumYears<-apply(A[,grepl("wage_sum_(1[23456])_Uncondit",names(A))],1,function(x) sum(x>0))
  A$Total_Wages30s_logUncond<-log(1+apply(A[,grepl("wage_sum_(1[23456])_Uncondit",names(A))],1,function(x) sum(x))/14)
  A$Total_Wages30s_logCondit<-log(1+apply(A[,grepl("wage_sum_(1[23456])_Uncondit",names(A))],1,function(x) sum(x))/A$Total_Wages30s_NumYears)
  A$Total_Wages_Condit<-apply(A[,grepl("wage_sum_([6-9]|1[0-6])_Uncondit",names(A))],1,function(x) sum(x))/A$Total_Wages_NumYears
  A$Total_Wages30s_Condit<-apply(A[,grepl("wage_sum_(1[23456])_Uncondit",names(A))],1,function(x) sum(x))/A$Total_Wages30s_NumYears
  A$Total_Wages30s_Num75000<-apply(A[,grepl("wage_sum_(1[23456])_Uncondit",names(A))],1,function(x) sum(x>75000))
  A$Total_Wages30s_Num100000<-apply(A[,grepl("wage_sum_(1[23456])_Uncondit",names(A))],1,function(x) sum(x>100000))
  A$Total_Wages30s_Num150000<-apply(A[,grepl("wage_sum_(1[23456])_Uncondit",names(A))],1,function(x) sum(x>150000))
  
  #Define a UC-only STEM variable
  stem<-unlist(read.csv("Data/Raw/Definition of STEM.csv",stringsAsFactors = F)) ; stem<-gsub("^.*[^0-9]([0-9]{6})[^0-9].*$","\\1",gsub("[.]","",stem[grepl("[.]",stem)])) ; stem[nchar(stem)==5]<-paste0("0",stem[nchar(stem)==5]) ;  for(z in stem[grepl("X",stem)]) stem<-c(stem,paste0(gsub("X","",z),addzeros(1:9999,4)))
    A$STEM_NoNSC<-NA_to_F(A$CIP6RECENT%in%stem) ; A$STEM_NoNSC[is.na(A$HASGRAD6)]<-NA
    A$STEM_NoNSC_Condit<-A$STEM_NoNSC ; A$STEM_NoNSC_Condit[NA_to_T(A$HASGRAD6==0)]<-NA

  A$AIS_UCSD<-round((A$hsgpa_censor*1000+A$SATIV+A$SATIM+A$SATII_W+A$SATII_M)/10)*10 #San Diego didn't include 3rd SATII in its AIS calculation
  for(s in c("01","03","04","05","06","07","08","09","10")){
    A[NA_to_F(A[,paste0("CMAJPRP",s)]=="149"),paste0("Major",s)]<-"Bioengineering"
  }
  A$AtoG<-A$HSCRSE_TYP_TOTSEMS_A>=4 & A$HSCRSE_TYP_TOTSEMS_B>=8 & A$HSCRSE_TYP_TOTSEMS_C>=6 & A$HSCRSE_TYP_TOTSEMS_D>=4 & A$HSCRSE_TYP_TOTSEMS_E>=4 & A$HSCRSE_TYP_TOTSEMS_F>=1 #Indicator for satisfying UC's A-to-F minimum high school course requirement
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces Figures I and B-1 to B-8.
  ##########################################################################################################
  ##########################################################################################################
  annuals<-T ; overall<-T
  A$Count<-1
  campuses<-c("BK","DV","LA","RV","SD","SC","SB","IR")
  count<-1 ; for(s in c("01","03","04","05","06","07","08","09")){
    A$ENR<-A[,paste0("ENR",s)]
    A$Major00<-paste(A[,paste0("Major",s)],A[,paste0("INTEND_FIELD_",campuses[count],"_STUDY")])
    A$ADM<-NA
    A$ADM[grepl("Not Retained|^Eligible - Other|Not Eligible - UC eligible, Deferred by College",A[,paste0("ELGBLITY_CMP_",campuses[count],"_DECN_DESC")])]<-0 #Take-up of 'other' appears tiny
    A$ADM[grepl("^Eligible - Special|^Eligible - REGULAR|^Eligible - CE",A[,paste0("ELGBLITY_CMP_",campuses[count],"_DECN_DESC")])]<-1
    #if(campuses[count]=="RV") A$ADM[grepl("No Action, retain, in progress",A[,paste0("ELGBLITY_CMP_",campuses[count],"_DECN_DESC")])]<-0 #Not sure if should do this; Riverside's basically just non-selective among eligible students
    A$APP<-!is.na(A$ADM)
    A$LS<-A[,paste0("LS",s)]|A[,paste0("CMAJPRP",s)]%in%c("","000","0")
    A$ais<-A$AIS #; if(s=="06") A$ais<-A$AIS_UCSD #Note: Not differentiating for now.
    aisrange<-4900:7900 #; if(s=="06") aisrange<-4400:7150
    if(annuals){
      for(i in 1994:2001){
        nu<-A[A$APP==1&A$adm_type=="Freshman"&A$LS&A$ais%in%aisrange&A$CETHNICA%in%c("P","E","D","F","G","H","N","V")&A$RESIDENCY_STATUS_DESC=="California Resident Student"&A$NotIneligible&A$YEARAPAY==i,]
        u<-A[A$APP==1&A$adm_type=="Freshman"&A$LS&A$ais%in%aisrange&A$CETHNICA%in%c("A","B","C")&A$RESIDENCY_STATUS_DESC=="California Resident Student"&A$NotIneligible&A$YEARAPAY==i,]
        hi<-c()
        for(z in aisrange[aisrange%%10==0]){
          hi<-rbind(hi,data.frame(
            ais=z,
            ADM_URM=mean(u$ADM[abs(u$ais-z)<=50])*100,
            ADM_NonURM=mean(nu$ADM[abs(nu$ais-z)<=50])*100,
            Count_URM=sum(abs(u$ais-z)<=50)/11,
            Count_NonURM=sum(abs(nu$ais-z)<=50)/11,
            ADM_URM_raw=mean(u$ADM[abs(u$ais-z)<=0])*100,
            ADM_NonURM_raw=mean(nu$ADM[abs(nu$ais-z)<=0])*100
          ))
        }
        tiff(paste0("Figures/HistAdm_",s,"_",i,".tiff"),6,5,"in",res=300)
          par(mar = c(5,4,1,4))
          plot(hi$ais,hi$ADM_NonURM,xlab="Academic Index",ylab="",ylim=c(0,100),pch=20,cex=.4,col="#0000FF",cex.axis=1.2,cex.lab=1.2)
          mtext(side = 2, line = 2.6, 'Percent Admitted (dots)',cex=1.2)
          lines(hi$ais,hi$Count_NonURM/nrow(nu)*5000,col="#9C9CFF") 
          lines(hi$ais,hi$Count_URM/nrow(u)*5000,col="#FF9E9E") 
          points(hi$ais,hi$ADM_URM,pch=18,cex=.4,col="#FF0000")
          points(hi$ais,hi$ADM_NonURM,pch=20,cex=.4,col="#0000FF")
          axis(side = 4,at=c(0,25,50,75,100),labels=c("0","0.5","1","1.5","2"),cex.axis=1.2)
          mtext(side = 4, line = 2.6, 'Percent of Applicants (lines)',cex=1.2)
        dev.off()
        tiff(paste0("Figures/HistAdm_",s,"_",i,"_NoDist.tiff"),5,5,"in",res=300)
          plot(hi$ais,hi$ADM_NonURM,xlab="Academic Index",ylab="",ylim=c(0,100),pch=20,cex=.4,col="#0000FF",cex.axis=1.2,cex.lab=1.2)
          points(hi$ais,hi$ADM_URM,pch=18,cex=.4,col="#FF0000")
          points(hi$ais,hi$ADM_NonURM,pch=20,cex=.4,col="#0000FF")
        dev.off()
        tiff(paste0("Figures/HistAdm_",s,"_",i,"_NoSmooth.tiff"),5,5,"in",res=300)
          plot(hi$ais,hi$ADM_NonURM,xlab="Academic Index",ylab="",ylim=c(0,100),pch=20,cex=.4,col="#0000FF",cex.axis=1.2,cex.lab=1.2)
          mtext(side = 2, line = 2.6, 'Percent Admitted (dots)',cex=1.2)
          lines(hi$ais,hi$Count_NonURM/nrow(nu)*5000,col="#9C9CFF") 
          lines(hi$ais,hi$Count_URM/nrow(u)*5000,col="#FF9E9E") 
          points(hi$ais,hi$ADM_URM_raw,pch=18,cex=.4,col="#FF0000")
          points(hi$ais,hi$ADM_NonURM_raw,pch=20,cex=.4,col="#0000FF")
          axis(side = 4,at=c(0,25,50,75,100),labels=c("0","0.5","1","1.5","2"),cex.axis=1.2)
          mtext(side = 4, line = 2.6, 'Percent of Applicants (lines)',cex=1.2)
        dev.off()
      }
    }
    if(overall){
      a<-data.frame()
      temp<-A[A$APP==1&A$adm_type=="Freshman"&A$LS&A$RESIDENCY_STATUS_DESC=="California Resident Student"&A$NotIneligible,]
      temp$YEARAPAY[temp$YEARAPAY%%2==1]<-temp$YEARAPAY[temp$YEARAPAY%%2==1]-1
      for(j in -5:5){ #Moving average, 50 points to either side
        temp$temp<-temp$ais+(10*j)
        a<-rbind(a,temp)
      } ; a$ais<-a$temp
      hi<-aggregate(ADM~ais+YEARAPAY,a[a$CETHNICA%in%c("P","E","D","F","G","H","N","V"),],mean)
      hi1<-aggregate(ADM~ais+YEARAPAY,a[a$CETHNICA%in%c("A","B","C","J"),],mean) #NOTE: Including J here, for completeness
      hi2<-aggregate(Count~ais+YEARAPAY,a[a$CETHNICA%in%c("A","B","C","J"),],sum)
      names(hi1)[3]<-"ADM_URM"
      hi<-merge(hi,hi1,all=T)
      hi<-merge(hi,hi2,all=T)
      hi$Count<-hi$Count/22
      hi$Diff<-(hi$ADM_URM-hi$ADM)*100
      hi_all<-hi ; hi<-hi[hi$ais%in%aisrange,]
      tiff(paste0("Figures/URMDiff_",s,".tiff"),5,3.5,"in",res=300)
        plot(hi$ais[hi$YEARAPAY==1994],hi$Diff[hi$YEARAPAY==1994],pch="",xlim=c(5000,7900),ylim=c(0,100),xlab="Academic Index",ylab="Percentage Points",yaxp=c(0,100,4),mgp=c(2.2,1,0))
        lines(hi$ais[hi$YEARAPAY==1994],hi$Diff[hi$YEARAPAY==1994],col="black",lwd=2)
        lines(hi$ais[hi$YEARAPAY==1996],hi$Diff[hi$YEARAPAY==1996],col="gray25")
        lines(hi$ais[hi$YEARAPAY==1998],hi$Diff[hi$YEARAPAY==1998],col="gray50",lty=2)
        lines(hi$ais[hi$YEARAPAY==2000],hi$Diff[hi$YEARAPAY==2000],col="gray75",lty=3)
        text(x=6850-700*(s=="06"&(0==1))-2000*(s=="01")+200*(!s%in%c("01","04")),y=100-4,pos = 4,labels=paste0("'94-5: ",ex(sum(hi_all$Diff[hi_all$YEARAPAY==1994]*hi_all$Count[hi_all$YEARAPAY==1994],na.rm=T)/100,0,justnum = T,comma=T)),cex = 1.4,font=2)
        text(x=6850-700*(s=="06"&(0==1))-2000*(s=="01")+200*(!s%in%c("01","04")),y=92-9,pos = 4,labels=paste0("'96-7: ",ex(sum(hi_all$Diff[hi_all$YEARAPAY==1996]*hi_all$Count[hi_all$YEARAPAY==1996],na.rm=T)/100,0,justnum = T,comma=T)),cex = 1.4,col="gray25")
        text(x=6850-700*(s=="06"&(0==1))-2000*(s=="01")+200*(!s%in%c("01","04")),y=84-14,pos = 4,labels=paste0("'98-9: ",ex(sum(hi_all$Diff[hi_all$YEARAPAY==1998]*hi_all$Count[hi_all$YEARAPAY==1998],na.rm=T)/100,0,justnum = T,comma=T)),cex = 1.4,col="gray50") 
        text(x=6850-700*(s=="06"&(0==1))-2000*(s=="01")+200*(!s%in%c("01","04")),y=76-19,pos = 4,labels=paste0("'00-1: ",ex(sum(hi_all$Diff[hi_all$YEARAPAY==2000]*hi_all$Count[hi_all$YEARAPAY==2000],na.rm=T)/100,0,justnum = T,comma=T)),cex = 1.4,col="gray75")
      dev.off()
    }
    count<-count+1
    print(s) ; gc()
  }
  #Legend
  tiff(paste0("Figures/Legend_AnnualURM.tiff"),width=8,height=3,units="in",res=500)
    plot(100000,100000,bty="n",box=F,yaxt="n",xaxt="n",xlim=c(0,50),ylim=c(0,10),xlab="",ylab="")
    legend(1,8,legend=c("Non-URM Admit Rate (left axis)"," ","Non-URM Distribution (right axis)","URM Admit Rate (left axis)"," ","URM Distribution (right axis)"),lty=c(NA,NA,1,NA,NA,1),pch=c(20,NA,NA,18,NA,NA),col=c("#0000FF",NA,"#9C9CFF","#FF0000",NA,"#FF9E9E"),text.width = c(13.8,13.8,13.8),cex=.6,ncol = 2)
  dev.off()
  
  
  
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces Figures A-1 and A-5.
  ##########################################################################################################
  ##########################################################################################################
  Results<-data.frame(matrix(NA,8,64),stringsAsFactors = F)
  count<-1 ; s<-"01"
  cs<-c("BK","DV","LA","RV","SD","SC","SB","IR")
  for(s in c("01","03","04","05","06","07","08","09")){
    countt<-0 ; for(type in c("","_overall")){
      A$Count<-1
      A$ADM<-NA
      A$ADM[grepl("Not Retained|^Eligible - Other|Not Eligible - UC eligible, Deferred by College",A[,paste0("ELGBLITY_CMP_",cs[count],"_DECN_DESC")])]<-0
      A$ADM[grepl("Eligible - Special|Eligible - REGULAR|Eligible - CE",A[,paste0("ELGBLITY_CMP_",cs[count],"_DECN_DESC")])]<-1
      if(type=="_overall") A$ADM<-A[,paste0("ADM",s)]
      A$APP<-!is.na(A$ADM)
      A$LS<-A[,paste0("LS",s)]|A[,paste0("CMAJPRP",s)]%in%c("","000","0")
      A$ais<-A$AIS ; if(s=="06") A$ais<-A$AIS_UCSD
      aisrange<-4500:8000 ; if(s=="06") aisrange<-4000:7200
      if(type=="") cond<-A$APP==1&A$adm_type=="Freshman"&A$LS&A$ais%in%aisrange&A$RESIDENCY_STATUS_DESC=="California Resident Student"&A$NotIneligible else cond<-A$APP==1&A$adm_type=="Freshman"&A$RESIDENCY_STATUS_DESC=="California Resident Student"
      
      A$EthGroup<-A$URM ; A$EthGroup[A$CETHNICA=="J"]<-2
      temp<-A[cond,]
      hi<-aggregate(.~ais+YEARAPAY+EthGroup,A[cond,c("ais","YEARAPAY","EthGroup","Count","ADM")],sum) ; names(hi)[(length(hi)-1):length(hi)]<-c("SumTotal","Admitted") ; temp<-merge(temp,hi,all=T)
      A<-A[,names(A)!=paste0("L",s)] ; A<-merge(A,hi,all.x=T) ; A[,paste0("L",s)]<-(A$Admitted-NA_to_F(A$ADM))/(A$SumTotal-NA_to_F(A$ADM)) ; 
      A[is.infinite(A[,paste0("L",s)])&NA_to_F(A$AIS>6500),paste0("L",s)]<-1 ; A[is.infinite(A[,paste0("L",s)]),paste0("L",s)]<-0
      A<-A[,!names(A)%in%c("SumTotal","Admitted")] 
      hi<-aggregate(.~SAT+YEARAPAY+EthGroup,A[cond,c("SAT","YEARAPAY","EthGroup","Count","ADM")],sum) ; names(hi)[(length(hi)-1):length(hi)]<-c("SumTotal_SAT","Admitted_SAT") ; temp<-merge(temp,hi,all=T)
      hi<-aggregate(.~hsgpa+YEARAPAY+EthGroup,A[cond,c("hsgpa","YEARAPAY","EthGroup","Count","ADM")],sum) ; names(hi)[(length(hi)-1):length(hi)]<-c("SumTotal_GPA","Admitted_GPA") ; temp<-merge(temp,hi,all=T)
      hi<-aggregate(.~ais+YEARAPAY,A[cond,c("ais","YEARAPAY","Count","ADM")],sum) ; names(hi)[(length(hi)-1):length(hi)]<-c("SumTotal_NoEth","Admitted_NoEth") ; temp<-merge(temp,hi,all=T)
      
      names(Results)[count*8-7+countt]<-paste0(cs[count],type) 
      temp$Likelihood<-(temp$Admitted-temp$ADM)/(temp$SumTotal-temp$ADM) ; temp$Likelihood[is.infinite(temp$Likelihood)]<-NA
      for(y in 1994:2001){
        Results[y-1993,paste0(cs[count],type)]<-round(summary(lm(ADM~Likelihood,data=temp[temp$YEARAPAY==y,]))$r.squared,3)
      }
      count1<-6
      for(v in c("_NoEth","_SAT","_GPA")){
        names(Results)[count*8-count1+countt]<-paste0(cs[count],v,type)
        temp$Likelihood<-(temp[,paste0("Admitted",v)]-temp$ADM)/(temp[,paste0("SumTotal",v)]-temp$ADM) ; temp<-temp[!is.infinite(temp$Likelihood),]
        for(y in 1994:2001){
          Results[y-1993,paste0(cs[count],v,type)]<-round(summary(lm(ADM~Likelihood,data=temp[temp$YEARAPAY==y,]))$r.squared,3)
        }
        count1<-count1-1
      }
      countt<-countt+4
    }
    print(cs[count])
    count<-count+1
  }
  save(Results,file="Figures/AA_Maimonides_R2_byCampus.Rda")
  load("Figures/AA_Maimonides_R2_byCampus.Rda")
  cs<-c("BK","DV","LA","RV","SD","SC","SB","IR")
  count<-1 ; for(s in c("01","03","04","05","06","07","08","09")){
    tiff(paste0("Figures/Maimonides_",s,".tiff"),width=500,height=500)
      plot(1994:2001,Results[,cs[count]],ylim=c(0,1),xlab="Year",ylab=expression(paste(R^2," on Admission")),cex=2,cex.axis=2,cex.lab=1.5,lty=1,type="l",lwd=2,mgp=c(2.5,1,0))
      lines(c(1995.5,1995.5),c(-20,20),col="gray80",lty=2)
      lines(c(1997.5,1997.5),c(-20,20),col="gray80",lty=2)
      lines(1994:2001,Results[,cs[count]],lwd=2)
      lines(1994:2001,Results[,paste0(cs[count],"_NoEth")])
      lines(1994:2001,Results[,paste0(cs[count],"_SAT")],lty=2)
      lines(1994:2001,Results[,paste0(cs[count],"_GPA")],lty=3)
      legend(1998,.98,lty=c(1,1,2,3),lwd=c(2,1,1,1),legend=c("AI + Eth.","AI","SAT","HS GPA"),cex=1.5)
    dev.off()
    count<-count+1
  }
  #Plot gaps for all eight campuses
  for(type in c("","_overall")){
    tiff(paste0("Figures/Maimonides_Diffs",type,".tiff"),width=5,height=5,res=300,units = "in")
      plot(1994:2001,Results$BK-Results$BK_NoEth,ylim=c(0,.4),xlab="Year",ylab=expression(paste("Difference in ",R^2)),cex=2,cex.lab=1.1,cex.axis=1.2,mgp=c(2.4,1,0),lty=1,type="l",lwd=2,yaxt="n")
      box()
      axis(2, labels = c("0","0.1","0.2","0.3","0.4","0.5"), at = seq(0,0.5,0.1),cex.axis=1.3)
      lines(c(1995.5,1995.5),c(-20,20),col="gray80",lty=2)
      lines(c(1997.5,1997.5),c(-20,20),col="gray80",lty=2)
      lines(1994:2001,Results[,paste0("BK",type)]-Results[,paste0("BK_NoEth",type)],lwd=3,col="black")
      lines(1994:2001,Results[,paste0("LA",type)]-Results[,paste0("LA_NoEth",type)],lwd=3,col="gray65")
      lines(1994:2001,Results[,paste0("SD",type)]-Results[,paste0("SD_NoEth",type)],lwd=3,col="gray30")
      lines(1994:2001,Results[,paste0("DV",type)]-Results[,paste0("DV_NoEth",type)],lwd=2,col="black",lty=3)
      lines(1994:2001,Results[,paste0("IR",type)]-Results[,paste0("IR_NoEth",type)],lwd=2,col="gray50",lty=3)
      lines(1994:2001,Results[,paste0("SB",type)]-Results[,paste0("SB_NoEth",type)],lwd=1,col="black",lty=2)
      lines(1994:2001,Results[,paste0("SC",type)]-Results[,paste0("SC_NoEth",type)],lwd=1,col="gray65",lty=2)
      lines(1994:2001,Results[,paste0("RV",type)]-Results[,paste0("RV_NoEth",type)],lwd=1,col="gray30",lty=2)
      legend(1998.25,.4,col=c("black","gray65","gray30","black","gray50","black","gray50"),lwd=c(3,3,3,2,2,1,1),lty=c(1,1,1,3,3,2,2),legend=c("Berkeley","UCLA","San Diego","Davis","Irvine","Santa Barbara","Santa Cruz"),cex=.8)
    dev.off()
  }
  
  
  
  #Test the UC Berkeley 96-97 non-URM admissions discontinuity
  A$ADM<-NA
    A$ADM[grepl("Not Retained|^Eligible - Other|Not Eligible - UC eligible, Deferred by College",A[,paste0("ELGBLITY_CMP_","BK","_DECN_DESC")])]<-0 #Take-up of 'other' appears tiny
    A$ADM[grepl("Eligible - Special|Eligible - REGULAR|Eligible - CE",A[,paste0("ELGBLITY_CMP_","BK","_DECN_DESC")])]<-1
  D<-A[!is.na(A$ADM)&A$LS01&A$YEARAPAY%in%c(1996,1997)&A$adm_type=="Freshman"&A$NotIneligible,]
    D$AIS[D$YEARAPAY==1996]<-D$AIS[D$YEARAPAY==1996]+70
  c<-7425
  
  #McCrary Test
  DCdensity(D$AIS[D$RESIDENCY_STATUS_DESC=="California Resident Student"&((D$CETHNICA%in%c("P","E","D","F","G","H","N","V")))],c=c,ext.out=T)[1:4] #Looks good
  hist(D$AIS[D$AIS%in%6010:8000],breaks = seq(6000,8000,10),main="",xlab="Academic Index")
    abline(v=7430,col="blue")
  
  cond<-D$CETHNICA%in%c("P","E","D","F","G","H","N","V","K","L","M")&D$AIS%in%seq(c-400,c+400)&D$AIS<7951
  condreg<-D$CETHNICA%in%c("P","E","D","F","G","H","N","V","K","L","M") ; D$Y97<-D$YEARAPAY==1997
  D$Above<-D$AIS>c
  D$AIS_Above<-(D$AIS-c)*D$Above ; D$AIS_Below<-(D$AIS-c)*(1-D$Above)
  D$AIS_Above2<-(D$AIS-c)^2*D$Above ; D$AIS_Below2<-(D$AIS-c)^2*(1-D$Above)
  D$Count<-1
  
  #Check that observables are balanced across the threshold
  hi<-lm(Total_Wages30s_logCondit~Female*factor(CETHNICA)+factor(educ_lvl_max)+logincome_parent_zeros+income_parent_zeros+factor(YEARAPAY),data=A[A$YEARAPAY%in%c(1996,1997)&A$adm_type=="Freshman"&A$NotIneligible&!A$AL_KEY%in%D$AL_KEY[cond],])
  D$predict_wage[!is.na(D$Total_Wages30s_logCondit)]<-predict(hi,D[!is.na(D$Total_Wages30s_logCondit),])
  reg<-rdrobust((D$predict_wage[cond]),D$AIS[cond],c=c,covs=D[cond,c("YEARAPAY")]) ; reg<-c(reg$coef[1],reg$se[1]) ; reg
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces Figure IX.
  ##########################################################################################################
  ##########################################################################################################
  #Panel a: Berkeley admission
  reg<-rdrobust(D$ADM01[condreg]*100,D$AIS[condreg],c=c,covs=D[condreg,c("YEARAPAY")]) ; reg<-c(reg$coef[1],reg$se[1])
  g<-rdplot(y=D$ADM01[cond]*100,x=D$AIS[cond],c=c,nbins=1000,p=1,kernel="tri",y.lim = c(0,100),col.lines="deepskyblue",x.label="Academic Index",y.label="Percent",x.lim=c(7030,7800),title="",col.dots = "black")[[2]]+ggplot2::annotate("text",x=7625,y=5,label=paste0("β = ",ex(reg[1],1,justnum=T,comma=T)," (",ex(reg[2],1,justnum=T,comma=T),")"),size=6)+
    theme(axis.text.y   = element_text(size=17.5),
          axis.text.x   = element_text(size=17.5),
          axis.title.y  = element_text(size=17.5),
          axis.title.x  = element_text(size=17.5),
          panel.background = element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),
          panel.border = element_rect(colour = "black", fill=NA, size=1)
    )
  ggsave(filename="BklyPre209_RD_ADM01.tiff",g,"tiff",path = "Figures",width = 5,height=3.5,units = "in",dpi = 300)
  
  #Panel b: Berkeley enrollment
  reg<-rdrobust(D$ENR01[condreg]*100,D$AIS[condreg],c=c,covs=D[condreg,c("YEARAPAY")]) ; reg<-c(reg$coef[1],reg$se[1])
  g<-rdplot(y=D$ENR01[cond]*100,x=D$AIS[cond],c=c,nbins=20,p=1,kernel="tri",y.lim = c(10,45),col.lines="deepskyblue",x.label="Academic Index",y.label="Percent",x.lim=c(7030,7800),title="",col.dots = "black")[[2]]+ggplot2::annotate("text",x=7625,y=11.75,label=paste0("β = ",ex(reg[1],1,justnum=T,comma=T)," (",ex(reg[2],1,justnum=T,comma=T),")"),size=6)+
    theme(axis.text.y   = element_text(size=17.5),
          axis.text.x   = element_text(size=17.5),
          axis.title.y  = element_text(size=17.5),
          axis.title.x  = element_text(size=17.5),
          panel.background = element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),
          panel.border = element_rect(colour = "black", fill=NA, size=1)
    )
  ggsave(filename="BklyPre209_RD_ENR01.tiff",g,"tiff",path = "Figures",width = 5,height=3.5,units = "in",dpi = 300)
  
  #Panel c: CFSTY wage value-added
  reg<-rdrobust(D$FE_WageVA_Chetty_AA[condreg]/1000,D$AIS[condreg],c=c,covs=D[condreg,c("YEARAPAY")]) ; reg<-c(reg$coef[1],reg$se[1])
  cond_spec<-D$CETHNICA%in%c("P","E","D","F","G","H","N","V","K","L","M")&D$AIS%in%seq(c-500,c+500)&D$AIS<7951 #Otherwise 'the leading minor of order 5 is not positive definite'. Doesn't change estimate.
  g<-rdplot(y=(D$FE_WageVA_Chetty_AA[cond_spec]-D$FE_WageVA_Chetty_AA[grepl("UNIVERSITY - LONG BEACH",D$NSC_EnrFirstIncCC_College_Name)][1])/1000,x=D$AIS[cond_spec],c=c,nbins=20,h=400,p=1,kernel="tri",y.lim = c(11,25),col.lines="deepskyblue",x.label="Academic Index",y.label="Dollars ('000s)",x.lim=c(7030,7800),title="",col.dots = "black")[[2]]+ggplot2::annotate("text",x=7625,y=11.7,label=paste0("β = ",ex(reg[1],2,justnum=T,comma=T)," (",ex(reg[2],2,justnum=T,comma=T),")"),size=6)+
    theme(axis.text.y   = element_text(size=17.5),
          axis.text.x   = element_text(size=17.5),
          axis.title.y  = element_text(size=17.5),
          axis.title.x  = element_text(size=17.5),
          panel.background = element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),
          panel.border = element_rect(colour = "black", fill=NA, size=1)
    )
  ggsave(filename="BklyPre209_RD_FE_WageVA_Chetty_AA.tiff",g,"tiff",path = "Figures",width = 5,height=3.5,units = "in",dpi = 300)  
  
  #Panel d: Grad degree
  reg<-rdrobust(D$MA[condreg]*100,D$AIS[condreg],c=c,covs=D[condreg,c("YEARAPAY")]) ; reg<-c(reg$coef[1],reg$se[1])
  g<-rdplot(y=(D$MA[cond]*100),x=D$AIS[cond],c=c,nbins=16,p=1,kernel="tri",y.lim = c(40,70),col.lines="deepskyblue",x.label="Academic Index",y.label="Percent",x.lim=c(7030,7800),title="",col.dots = "black")[[2]]+ggplot2::annotate("text",x=7625,y=41.5,label=paste0("β = ",ex(reg[1],1,justnum=T,comma=T)," (",ex(reg[2],1,justnum=T,comma=T),")"),size=6)+
    theme(axis.text.y   = element_text(size=17.5),
          axis.text.x   = element_text(size=17.5),
          axis.title.y  = element_text(size=17.5),
          axis.title.x  = element_text(size=17.5),
          panel.background = element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),
          panel.border = element_rect(colour = "black", fill=NA, size=1)
    )
  ggsave(filename="BklyPre209_RD_MA.tiff",g,"tiff",path = "Figures",width = 5,height=3.5,units = "in",dpi = 300)
  
  #Panel e: Average log early-30s wages
  reg<-rdrobust(D$Total_Wages30s_logCondit[condreg],D$AIS[condreg],c=c,covs=D[condreg,c("YEARAPAY")]) ; reg<-c(reg$coef[1],reg$se[1],reg$bws[1,1])
  cond_spec<-D$CETHNICA%in%c("P","E","D","F","G","H","N","V","K","L","M")&D$AIS%in%seq(c-600,c+600)&D$AIS<7951 #Otherwise 'the leading minor of order 5 is not positive definite'. Doesn't change estimate.
  g<-rdplot(y=(D$Total_Wages30s_logCondit[cond_spec]),x=D$AIS[cond_spec],c=c,nbins=24,p=1,kernel="triangular",h=400,y.lim = c(10.5,11.5),col.lines="deepskyblue",x.label="Academic Index",y.label="Log Dollars",x.lim=c(7030,7800),title="",col.dots = "black")[[2]]+ggplot2::annotate("text",x=7625,y=10.55,label=paste0("β = ",ex(reg[1],2,justnum=T,comma=T)," (",ex(reg[2],2,justnum=T,comma=T),")"),size=6)+
    theme(axis.text.y   = element_text(size=17.5),
          axis.text.x   = element_text(size=17.5),
          axis.title.y  = element_text(size=17.5),
          axis.title.x  = element_text(size=17.5),
          panel.background = element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),
          panel.border = element_rect(colour = "black", fill=NA, size=1)
    )
  ggsave(filename="BklyPre209_RD_Total_Wages30s_logCondit.tiff",g,"tiff",path = "Figures",width = 5,height=3.5,units = "in",dpi = 300) 
  
  #Panel f: Number of $150k+ years in early 30s
  reg<-rdrobust(D$Total_Wages30s_Num150000[condreg],D$AIS[condreg],c=c,covs=D[condreg,c("YEARAPAY")]) ; reg<-c(reg$coef[1],reg$se[1])
  g<-rdplot(y=(D$Total_Wages30s_Num150000[cond]),x=D$AIS[cond],c=c,nbins=16,p=1,kernel="tri",y.lim = c(0.35,1.07),col.lines="deepskyblue",x.label="Academic Index",y.label="# of Years",x.lim=c(7030,7800),title="",col.dots = "black")[[2]]+ggplot2::annotate("text",x=7625,y=0.386,label=paste0("β = ",ex(reg[1],2,justnum=T,comma=T)," (",ex(reg[2],2,justnum=T,comma=T),")"),size=6)+
    theme(axis.text.y   = element_text(size=17.5),
          axis.text.x   = element_text(size=17.5),
          axis.title.y  = element_text(size=17.5),
          axis.title.x  = element_text(size=17.5),
          panel.background = element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),
          panel.border = element_rect(colour = "black", fill=NA, size=1)
    )
  ggsave(filename="BklyPre209_RD_Total_Wages30s_Num150000.tiff",g,"tiff",path = "Figures",width = 5,height=3.5,units = "in",dpi = 300) 
  
  
  
  
  
  ###
  #  UC Davis Discontinuity
  ###
  A$ADM<-NA
    A$ADM[grepl("Not Retained|^Eligible - Other|Not Eligible - UC eligible, Deferred by College",A[,paste0("ELGBLITY_CMP_","DV","_DECN_DESC")])]<-0 #Take-up of 'other' appears tiny
    A$ADM[grepl("Eligible - Special|Eligible - REGULAR|Eligible - CE",A[,paste0("ELGBLITY_CMP_","DV","_DECN_DESC")])]<-1
  D<-A[!is.na(A$ADM)&A$LS03&A$YEARAPAY%in%c(1996,1997)&A$adm_type=="Freshman"&A$NotIneligible,]
  D$AIS[D$YEARAPAY==1996]<-D$AIS[D$YEARAPAY==1996]+250
  c<-6245
  
  #McCrary Test
  DCdensity(D$AIS[D$RESIDENCY_STATUS_DESC=="California Resident Student"&((D$CETHNICA%in%c("P","E","D","F","G","H","N","V")))],c=c,ext.out=T)[1:4] #Fails: p=0.014
    DCdensity(D$AIS[D$RESIDENCY_STATUS_DESC=="California Resident Student"&((D$CETHNICA%in%c("P","E","D","F","G","H","N","V")))&D$YEARAPAY==1996],c=c,ext.out=T)[1:4] #Fails: p=0.016
    DCdensity(D$AIS[D$RESIDENCY_STATUS_DESC=="California Resident Student"&((D$CETHNICA%in%c("P","E","D","F","G","H","N","V")))&D$YEARAPAY==1997],c=c,ext.out=T)[1:4] #Fails: p=0.025
  tiff("Figures/DavisPre209_McCrary.tiff",5,5,"in",res=300)
    DCdensity(D$AIS[D$RESIDENCY_STATUS_DESC=="California Resident Student"&((D$CETHNICA%in%c("P","E","D","F","G","H","N","V")))],c=c,ext.out=T,bin = 30)[1:4]
    abline(v=6250,col="blue")
  dev.off()
  
  cond<-D$CETHNICA%in%c("P","E","D","F","G","H","N","V","K","L","M")&D$AIS%in%seq(c-400,c+400)&D$AIS<7951
  condreg<-D$CETHNICA%in%c("P","E","D","F","G","H","N","V","K","L","M") ; D$Y97<-D$YEARAPAY==1997
  D$Above<-D$AIS>c
  D$AIS_Above<-(D$AIS-c)*D$Above ; D$AIS_Below<-(D$AIS-c)*(1-D$Above)
  D$AIS_Above2<-(D$AIS-c)^2*D$Above ; D$AIS_Below2<-(D$AIS-c)^2*(1-D$Above)
  D$Count<-1
  
  #Check that observables are balanced across the threshold
  hi<-lm(Total_Wages30s_logCondit~Female*factor(CETHNICA)+factor(educ_lvl_max)+logincome_parent_zeros+income_parent_zeros+factor(YEARAPAY),data=A[A$YEARAPAY%in%c(1996,1997)&A$adm_type=="Freshman"&A$NotIneligible&!A$AL_KEY%in%D$AL_KEY[cond],])
  D$predict_wage[!is.na(D$Total_Wages30s_logCondit)]<-predict(hi,D[!is.na(D$Total_Wages30s_logCondit),])
  reg<-rdrobust((D$predict_wage[cond]),D$AIS[cond],c=c,covs=D[cond,c("YEARAPAY")]) ; reg<-c(reg$coef[1],reg$se[1]) ; reg
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces Figure J-1.
  ##########################################################################################################
  ##########################################################################################################
  #Panel a: Davis admission
  reg<-rdrobust(D$ADM03[condreg]*100,D$AIS[condreg],c=c,covs=D[condreg,c("YEARAPAY")]) ; reg<-c(reg$coef[1],reg$se[1])
  g<-rdplot(y=D$ADM03[cond]*100,x=D$AIS[cond],c=c,nbins=1000,p=1,kernel="tri",y.lim = c(0,100),col.lines="deepskyblue",x.label="Academic Index",y.label="Percent",x.lim=c(5850,6640),title="",col.dots = "black")[[2]]+ggplot2::annotate("text",x=6460,y=5,label=paste0("β = ",ex(reg[1],1,justnum=T,comma=T)," (",ex(reg[2],1,justnum=T,comma=T),")"),size=6)+
    theme(axis.text.y   = element_text(size=17.5),
          axis.text.x   = element_text(size=17.5),
          axis.title.y  = element_text(size=17.5),
          axis.title.x  = element_text(size=17.5),
          panel.background = element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),
          panel.border = element_rect(colour = "black", fill=NA, size=1)
    )
  ggsave(filename="DavisPre209_RD_ADM03.tiff",g,"tiff",path = "Figures",width = 5,height=3.5,units = "in",dpi = 300)
  
  #Panel b: Davis enrollment
  reg<-rdrobust(D$ENR03[condreg]*100,D$AIS[condreg],c=c,covs=D[condreg,c("YEARAPAY")]) ; reg<-c(reg$coef[1],reg$se[1])
  g<-rdplot(y=D$ENR03[cond]*100,x=D$AIS[cond],c=c,nbins=20,p=1,kernel="tri",y.lim = c(0,45),col.lines="deepskyblue",x.label="Academic Index",y.label="Percent",x.lim=c(5850,6640),title="",col.dots = "black")[[2]]+ggplot2::annotate("text",x=6460,y=2.25,label=paste0("β = ",ex(reg[1],1,justnum=T,comma=T)," (",ex(reg[2],1,justnum=T,comma=T),")"),size=6)+
    theme(axis.text.y   = element_text(size=17.5),
          axis.text.x   = element_text(size=17.5),
          axis.title.y  = element_text(size=17.5),
          axis.title.x  = element_text(size=17.5),
          panel.background = element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),
          panel.border = element_rect(colour = "black", fill=NA, size=1)
    )
  ggsave(filename="DavisPre209_RD_ENR03.tiff",g,"tiff",path = "Figures",width = 5,height=3.5,units = "in",dpi = 300)
  
  #Panel c: CFSTY wage value-added
  reg<-rdrobust(D$FE_WageVA_Chetty_AA[condreg]/1000,D$AIS[condreg],c=c,covs=D[condreg,c("YEARAPAY")]) ; reg<-c(reg$coef[1],reg$se[1])
  cond_spec<-D$CETHNICA%in%c("P","E","D","F","G","H","N","V","K","L","M")&D$AIS%in%seq(c-500,c+500)&D$AIS<7951 #Otherwise 'the leading minor of order 5 is not positive definite'. Doesn't change estimate.
  g<-rdplot(y=(D$FE_WageVA_Chetty_AA[cond_spec]-D$FE_WageVA_Chetty_AA[grepl("UNIVERSITY - LONG BEACH",D$NSC_EnrFirstIncCC_College_Name)][1])/1000,x=D$AIS[cond_spec],c=c,nbins=20,h=400,p=1,kernel="tri",y.lim = c(5,11),col.lines="deepskyblue",x.label="Academic Index",y.label="Dollars ('000s)",x.lim=c(5850,6640),title="",col.dots = "black")[[2]]+ggplot2::annotate("text",x=6460,y=5+(11-5)*0.05,label=paste0("β = ",ex(reg[1],2,justnum=T,comma=T)," (",ex(reg[2],2,justnum=T,comma=T),")"),size=6)+
    theme(axis.text.y   = element_text(size=17.5),
          axis.text.x   = element_text(size=17.5),
          axis.title.y  = element_text(size=17.5),
          axis.title.x  = element_text(size=17.5),
          panel.background = element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),
          panel.border = element_rect(colour = "black", fill=NA, size=1)
    )
  ggsave(filename="DavisPre209_RD_FE_WageVA_Chetty_AA.tiff",g,"tiff",path = "Figures",width = 5,height=3.5,units = "in",dpi = 300)  
  
  #Panel d: Grad degree
  reg<-rdrobust(D$MA[condreg]*100,D$AIS[condreg],c=c,covs=D[condreg,c("YEARAPAY")]) ; reg<-c(reg$coef[1],reg$se[1])
  g<-rdplot(y=(D$MA[cond]*100),x=D$AIS[cond],c=c,nbins=16,p=1,kernel="tri",y.lim = c(15,45),col.lines="deepskyblue",x.label="Academic Index",y.label="Percent",x.lim=c(5850,6640),title="",col.dots = "black")[[2]]+ggplot2::annotate("text",x=6460,y=15+30*0.05,label=paste0("β = ",ex(reg[1],1,justnum=T,comma=T)," (",ex(reg[2],1,justnum=T,comma=T),")"),size=6)+
    theme(axis.text.y   = element_text(size=17.5),
          axis.text.x   = element_text(size=17.5),
          axis.title.y  = element_text(size=17.5),
          axis.title.x  = element_text(size=17.5),
          panel.background = element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),
          panel.border = element_rect(colour = "black", fill=NA, size=1)
    )
  ggsave(filename="DavisPre209_RD_MA.tiff",g,"tiff",path = "Figures",width = 5,height=3.5,units = "in",dpi = 300)
  
  #Panel e: Average log early-30s wages
  reg<-rdrobust(D$Total_Wages30s_logCondit[condreg],D$AIS[condreg],c=c,covs=D[condreg,c("YEARAPAY")]) ; reg<-c(reg$coef[1],reg$se[1],reg$bws[1,1])
  cond_spec<-D$CETHNICA%in%c("P","E","D","F","G","H","N","V","K","L","M")&D$AIS%in%seq(c-600,c+600)&D$AIS<7951 #Otherwise 'the leading minor of order 5 is not positive definite'. Doesn't change estimate.
  g<-rdplot(y=(D$Total_Wages30s_logCondit[cond_spec]),x=D$AIS[cond_spec],c=c,nbins=24,p=1,kernel="triangular",h=400,y.lim = c(10.5,11.25),col.lines="deepskyblue",x.label="Academic Index",y.label="Log Dollars",x.lim=c(5850,6640),title="",col.dots = "black")[[2]]+ggplot2::annotate("text",x=6460,y=10.5+0.05*0.75,label=paste0("β = ",ex(reg[1],2,justnum=T,comma=T)," (",ex(reg[2],2,justnum=T,comma=T),")"),size=6)+
    theme(axis.text.y   = element_text(size=17.5),
          axis.text.x   = element_text(size=17.5),
          axis.title.y  = element_text(size=17.5),
          axis.title.x  = element_text(size=17.5),
          panel.background = element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),
          panel.border = element_rect(colour = "black", fill=NA, size=1)
    )
  ggsave(filename="DavisPre209_RD_Total_Wages30s_logCondit.tiff",g,"tiff",path = "Figures",width = 5,height=3.5,units = "in",dpi = 300) 
  
  #Panel f: Number of $150k+ years in early 30s
  reg<-rdrobust(D$Total_Wages30s_Num150000[condreg],D$AIS[condreg],c=c,covs=D[condreg,c("YEARAPAY")]) ; reg<-c(reg$coef[1],reg$se[1])
  g<-rdplot(y=(D$Total_Wages30s_Num150000[cond]),x=D$AIS[cond],c=c,nbins=16,p=1,kernel="tri",y.lim = c(0.1,.55),col.lines="deepskyblue",x.label="Academic Index",y.label="# of Years",x.lim=c(5850,6640),title="",col.dots = "black")[[2]]+ggplot2::annotate("text",x=6460,y=.1+0.05*0.45,label=paste0("β = ",ex(reg[1],2,justnum=T,comma=T)," (",ex(reg[2],2,justnum=T,comma=T),")"),size=6)+
    theme(axis.text.y   = element_text(size=17.5),
          axis.text.x   = element_text(size=17.5),
          axis.title.y  = element_text(size=17.5),
          axis.title.x  = element_text(size=17.5),
          panel.background = element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),
          panel.border = element_rect(colour = "black", fill=NA, size=1)
    )
  ggsave(filename="DavisPre209_RD_Total_Wages30s_Num150000.tiff",g,"tiff",path = "Figures",width = 5,height=3.5,units = "in",dpi = 300) 
  
  
}










if(AppChanges){
  load(paste0(secure_derived,"AffAct_Data.Rda"))
  a<-A[!is.na(A$HS_County)&A$CQUARTERAP==2&A$SCHTYPE=="A"&A$adm_type=="Freshman"&A$SEX%in%c("F","M"),]
  
  
  
  #Now build out from the school dataset
  G<-read.csv("Data/Derived/Graduation.csv",stringsAsFactors = F)[,-1]
    for(v in names(G)[grepl("Asian",names(G))]) G[,v]<-G[,v]-G[,gsub("Asian","Filipino",v)]
  
  G$HS_County[nchar(G$CDS_CODE)==13]<-as.integer(substr(G$CDS_CODE[nchar(G$CDS_CODE)==13],1,1))
    G$HS_County[nchar(G$CDS_CODE)==14]<-as.integer(substr(G$CDS_CODE[nchar(G$CDS_CODE)==14],1,2))
  G$HS_District[nchar(G$CDS_CODE)==13]<-as.integer(substr(G$CDS_CODE[nchar(G$CDS_CODE)==13],2,6))
    G$HS_District[nchar(G$CDS_CODE)==14]<-as.integer(substr(G$CDS_CODE[nchar(G$CDS_CODE)==14],3,7))
  G$HS_School[nchar(G$CDS_CODE)==13]<-as.integer(substr(G$CDS_CODE[nchar(G$CDS_CODE)==13],7,13))
    G$HS_School[nchar(G$CDS_CODE)==14]<-as.integer(substr(G$CDS_CODE[nchar(G$CDS_CODE)==14],8,14))
  G$YEARAPAY<-G$Year
  
  #Merge in number of apps
  a$AppAll<-a$NotIneligible ; a$AppSel<-(a$APP01==1|a$APP04==1)&a$NotIneligible
  a1<-aggregate(.~HS_County+HS_District+HS_School+YEARAPAY+EthGen,a[,c("HS_County","HS_District","HS_School","YEARAPAY","EthGen","AppAll","AppSel")],sum)
    a_1<-dcast(a1[,!names(a1)=="AppSel"],HS_County+HS_District+HS_School+YEARAPAY~EthGen) ; names(a_1)[5:ncol(a_1)]<-paste0(names(a_1)[5:ncol(a_1)],"_UCApp")
    a_2<-dcast(a1[,!names(a1)=="AppAll"],HS_County+HS_District+HS_School+YEARAPAY~EthGen) ; names(a_2)[5:ncol(a_2)]<-paste0(names(a_2)[5:ncol(a_2)],"_UCSelApp")
    a1<-merge(a_1,a_2)
  G<-merge(G,a1,all.x=T)
  for(v in names(a1)[grepl("_[MF]",names(a1))]) G[is.na(G[,v])&paste(G$HS_County,G$HS_District,G$HS_School)%in%unique(paste(a$HS_County,a$HS_District,a$HS_School)),v]<-0
  
  #Merge summed AIS score
  a$AIS_Count<-!is.na(a$AIS) ; a$AIS_Count[!a$NotIneligible]<-NA
  a1<-aggregate(.~HS_County+HS_District+HS_School+YEARAPAY+EthGen,a[a$NotIneligible,c("HS_County","HS_District","HS_School","YEARAPAY","EthGen","AIS","AIS_Count")],sum,na.rm=T)
    a_1<-dcast(a1[,!names(a1)=="AIS_Count"],HS_County+HS_District+HS_School+YEARAPAY~EthGen) ; names(a_1)[5:ncol(a_1)]<-paste0(names(a_1)[5:ncol(a_1)],"_AIS")
    a_2<-dcast(a1[,!names(a1)=="AIS"],HS_County+HS_District+HS_School+YEARAPAY~EthGen) ; names(a_2)[5:ncol(a_2)]<-paste0(names(a_2)[5:ncol(a_2)],"_AIS_Count")
    a1<-merge(a_1,a_2)
  G<-merge(G,a1,all.x=T)
  for(v in names(a1)[grepl("_[MF]",names(a1))]) G[is.na(G[,v])&paste(G$HS_County,G$HS_District,G$HS_School)%in%unique(paste(a$HS_County,a$HS_District,a$HS_School)),v]<-0 #Don't do this for SAT/HSGPA, but it works for summed AIS
  
  for(v in names(G)[grepl("_M",names(G))]) G[,gsub("_M","",v)]<-G[,v]+G[,gsub("_M","_F",v)] #Sum the male and female numbers
  
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces Figures VII and A-14.
  ##########################################################################################################
  ##########################################################################################################
  a$APP14<-(a$APP01==1|a$APP04==1)
  a$adm01<-NA
  a$adm01[grepl("Not Retained|^Eligible - Other|Not Eligible - UC eligible, Deferred by College",a[,paste0("ELGBLITY_CMP_","BK","_DECN_DESC")])]<-0 
  a$adm01[grepl("Eligible - Special|Eligible - REGULAR|Eligible - CE",a[,paste0("ELGBLITY_CMP_","BK","_DECN_DESC")])]<-1
  a$adm04<-NA
  a$adm04[grepl("Not Retained|^Eligible - Other|Not Eligible - UC eligible, Deferred by College",a[,paste0("ELGBLITY_CMP_","LA","_DECN_DESC")])]<-0 
  a$adm04[grepl("Eligible - Special|Eligible - REGULAR|Eligible - CE",a[,paste0("ELGBLITY_CMP_","LA","_DECN_DESC")])]<-1
  a$adm14[!is.na(a$adm01)|!is.na(a$adm04)]<-0 ; a$adm14[NA_to_F(a$adm01==1)|NA_to_F(a$adm04==1)]<-1
  campuses<-c("BK","DV","LA","RV","SD","SC","SB","IR")
  a$admAllUC<-NA
  for(j in 1:8) a$admAllUC[grepl("Not Retained|^Eligible - Other|Not Eligible - UC eligible, Deferred by College",a[,paste0("ELGBLITY_CMP_",campuses[j],"_DECN_DESC")])]<-0 #Take-up of 'other' appears tiny
  for(j in 1:8) a$admAllUC[grepl("Eligible - Special|Eligible - REGULAR|Eligible - CE",a[,paste0("ELGBLITY_CMP_",campuses[j],"_DECN_DESC")])]<-1
  a$APPAllUC<-1
  a$CETHNICA_Cat<-a$CETHNICA_Cat_Pred
  for(v in names(G)[grepl("_UC$",names(G))]) G[,gsub("_UC","_NotUC",v)]<-pmax(G[,gsub("_UC","",v)]-G[,v],0)
  
  for(elig in c("","_Inelig")){
    es<-c("Hispanic","Black","Asian")
    for(type in c("AllUC","01","04","14")){
      Ests<-data.frame(matrix(NA,15,13)) ; r<-1
      for(i in seq(4000,7800,200)){
        Ests[r,1]<-i
        a$elig<-a$NotIneligible ; if(elig=="_Inelig") a$elig<-!a$NotIneligible
        notuc<-"" ; if(elig=="_Inelig") notuc<-"Not"
        if(!is.na(as.integer(type))) temp<-a[a$AIS%in%(i:(i+199))&a$elig&a[,paste0("APP",type)]==1,] else temp<-a[a$AIS%in%(i:(i+199))&a$elig,]
        
        if(nrow(temp)<10){
          for(j in 1:3){
            Ests[r,j*3-1]<-0
            Ests[r,j*3]<-0.01
            Ests[r,j*3+1]<-1
          }
          r<-r+1 ; print(i)
          next
        }
        
        temp$Count<-1
        a1<-aggregate(.~HS_County+HS_District+HS_School+YEARAPAY+EthGen,temp[,c("HS_County","HS_District","HS_School","YEARAPAY","EthGen","Count")],sum,na.rm=T)
        a1_temp<-aggregate(.~HS_County+HS_District+HS_School+YEARAPAY+CETHNICA_Cat,temp[,c("HS_County","HS_District","HS_School","YEARAPAY","CETHNICA_Cat","Count")],sum,na.rm=T) ; names(a1_temp)[5]<-"EthGen" ; a1<-rbind(a1,a1_temp)
        a2<-melt(G[G$YEARAPAY%in%1994:2001,],id.vars=c("HS_County","HS_District","HS_School","YEARAPAY"),measure.vars=names(G)[grepl(paste0("_",notuc,"UC$"),names(G))])
        names(a2)[5:6]<-c("EthGen","Grads") ; a2$EthGen<-gsub(paste0("_",notuc,"UC"),"",a2$EthGen)
        a2<-a2[paste(a2$HS_County,a2$HS_District,a2$HS_School)%in%paste(a1$HS_County,a1$HS_District,a1$HS_School),]
        temp<-merge(a2,a1,all.x=T)
        temp$Count[is.na(temp$Count)]<-0 #304 cases of 0 grads & pos. UC grads. 319 where Grads<Count. Just drop...
        temp<-temp[temp$Grads!=0 & temp$Grads>=temp$Count,]
        temp$PercApply<-temp$Count/temp$Grads
        #For now, splitting by gender
        temp$FE1<-paste(temp$HS_County,temp$HS_District,temp$HS_School,temp$EthGen)
        temp$FE2<-paste(temp$HS_County,temp$HS_District,temp$HS_School,temp$YEARAPAY)
        temp$School<-paste(temp$HS_County,temp$HS_District,temp$HS_School)
        #for(y in c(1994,1996:2001)) for(e in c("Asian","Black","Hispanic","Filipino")) temp[,paste0(e,y)]<-temp$YEARAPAY==y&grepl(e,temp$EthGen)
        for(y in c(1994,1996,1998,2000)) for(e in c("Asian","Black","Hispanic")) temp[,paste0(e,y)]<-temp$YEARAPAY%in%c(y,y+1)&grepl(e,temp$EthGen)
        reg<-summary(felm(PercApply~Asian1996+Asian1998+Asian2000+
                            Hispanic1996+Hispanic1998+Hispanic2000+
                            Black1996+Black1998+Black2000|FE1+FE2+YEARAPAY|0|School,temp[grepl("(Asian|Black|Hispanic|White)$",temp$EthGen),],weights = temp$Grads[grepl("(Asian|Black|Hispanic|White)$",temp$EthGen)]),robust=T) #_[MF] or $
        for(j in 1:3){
          Ests[r,j*3-1]<-reg$coefficients[paste0(es[j],"1998TRUE"),1]*sum(G[G$Year%in%c(1998,1999),paste0(es[j],"_",notuc,"UC")])/2
          Ests[r,j*3]<-reg$coefficients[paste0(es[j],"1998TRUE"),2]*sum(G[G$Year%in%c(1998,1999),paste0(es[j],"_",notuc,"UC")])/2
          Ests[r,j*3+1]<-reg$coefficients[paste0(es[j],"1998TRUE"),4]
        }
        r<-r+1 ; print(i)
      }
      if(elig==""&type=="AllUC") save(Ests,file="Figures/AppDecline_AllUC_Ests.Rda")
      
      ys<-list(c(-300,150),c(-100,50),c(-300,300),c(-80,80))
      for(i in 1:3){ #4
        adms<-c() ; for(j in Ests[,1]) adms<-c(adms,mean(a[a$NotIneligible&a$CETHNICA_Cat==es[i]&a$AIS%in%j:(j+199)&a$YEARAPAY%in%c(1998,1999),paste0("adm",type)],na.rm=T))
        tiff(paste0("Figures/AppDecline_AIS_",type,"_",es[i],elig,".tiff"),width=1000,height=500)
          barplot(Ests[,i*3-1],
                  ylim=ys[[i]],
                  pch=19, xlab="Academic Index", ylab="Number of Students",
                  names.arg=Ests$X1,col="gray20",cex=2.5,cex.axis=1.8,cex.lab = 2,mgp=c(2.6,1,0)
          )
          barplot(Ests[,i*3-1]*adms, col='gray70', add=T,cex.axis=1.8,cex=2.5)
          arrows(seq(1,200,1.204)[1:20]-.295, Ests[,i*3-1]-Ests[,i*3]*1.96,seq(1,200,1.204)[1:20]-.295, Ests[,i*3-1]+Ests[,i*3]*1.96, length=0.05, angle=90, code=3,col="grey60",lty=2)
          text(20,ys[[i]][1]+(.15*(ys[[i]][2]-ys[[i]][1])),paste0("Apps: ",
                                                                  ex(sum(Ests[,i*3-1]),0,justnum = T,comma=T)," (",
                                                                  ex(sum(Ests[,i*3-1])/sum(a[a$CETHNICA_Cat==es[i]&a$YEARAPAY%in%c(1998,1999),paste0("APP",type)],na.rm=T)*100,1,justnum = T), #a$elig&
                                                                  "%)"),cex=2.5,pos = 3) 
          text(20,ys[[i]][1]+(.05*(ys[[i]][2]-ys[[i]][1])),paste0("Admits: ",
                                                                  ex(sum(Ests[,i*3-1]*NA_to_F(adms)),0,justnum = T,comma=T)," (",
                                                                  ex(sum(Ests[,i*3-1]*NA_to_F(adms))/sum(a[a$CETHNICA_Cat==es[i]&a$YEARAPAY%in%c(1998,1999),paste0("adm",type)],na.rm=T)*100,1,justnum = T), #a$elig&
                                                                  "%)"),
               cex=2.5,pos = 3)
          legend(0,ys[[i]][1]+(.3*(ys[[i]][2]-ys[[i]][1])),fill=c("black","gray70"),legend=c("Change in Apps","Apps Likely Admitted"),cex=2.2)
        dev.off()
      }
    }
  }
  
  
  
  
  
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces the statistics presented in Table A-15.
  ##########################################################################################################
  ##########################################################################################################
  R<-data.frame(matrix("",70,6),stringsAsFactors = F)
  for(type in c("AllUC","14")){
    if(!is.na(as.integer(type))) temp<-a[a$NotIneligible&a[,paste0("APP",type)]==1,] else temp<-a[a$NotIneligible,]
    temp$Count<-1
    a1<-aggregate(.~HS_County+HS_District+HS_School+YEARAPAY+EthGen,temp[,c("HS_County","HS_District","HS_School","YEARAPAY","EthGen","Count")],sum,na.rm=T)
    a1_temp<-aggregate(.~HS_County+HS_District+HS_School+YEARAPAY+CETHNICA_Cat,temp[,c("HS_County","HS_District","HS_School","YEARAPAY","CETHNICA_Cat","Count")],sum,na.rm=T) ; names(a1_temp)[5]<-"EthGen" ; a1<-rbind(a1,a1_temp)
    a2<-melt(G[G$YEARAPAY%in%1994:2001,],id.vars=c("HS_County","HS_District","HS_School","YEARAPAY"),measure.vars=names(G)[grepl("_UC$",names(G))])
    names(a2)[5:6]<-c("EthGen","Grads") ; a2$EthGen<-gsub("_UC","",a2$EthGen)
    a2<-a2[paste(a2$HS_County,a2$HS_District,a2$HS_School)%in%paste(a1$HS_County,a1$HS_District,a1$HS_School),]
    temp<-merge(a2,a1,all.x=T)
    temp$Count[is.na(temp$Count)]<-0 #304 cases of 0 grads & pos. UC grads. 319 where Grads<Count. Just drop...
    temp<-temp[temp$Grads!=0 & temp$Grads>=temp$Count,]
    temp$PercApply<-temp$Count/temp$Grads
    #For now, splitting by gender
    temp$FE1<-paste(temp$HS_County,temp$HS_District,temp$HS_School,temp$EthGen)
    temp$FE2<-paste(temp$HS_County,temp$HS_District,temp$HS_School,temp$YEARAPAY)
    temp$School<-paste(temp$HS_County,temp$HS_District,temp$HS_School)
    for(y in c(1995:2001)) for(e in c("Asian","Black","Hispanic","Filipino")) temp[,paste0(e,y)]<-temp$YEARAPAY==y&grepl(e,temp$EthGen)
    tf<-function(R,c,reg){
      r<-reg$coefficients[c(paste0("Black",c(1995:2001),"TRUE"),paste0("Hispanic",c(1995:2001),"TRUE"),paste0("Asian",c(1995:2001),"TRUE")),]
      R[seq(1,63,3),c]<-ex(r[,1],3)
      R[seq(2,64,3),c]<-ex(r[,2],3,par=T)
      
      R[66,c]<-"X"
      R[67,c]<-"X"
      
      R[69,c]<-ex(reg$r.squared,2)
      R[70,c]<-ex(reg$N,0,comma=T)
      return(R)
    }
    R<-tf(R,1+3*(type=="14"),summary(felm(PercApply~Asian1995+Asian1996+Asian1997+Asian1998+Asian1999+Asian2000+Asian2001+
                                            Black1995+Black1996+Black1997+Black1998+Black1999+Black2000+Black2001+
                                            Hispanic1995+Hispanic1996+Hispanic1997+Hispanic1998+Hispanic1999+Hispanic2000+Hispanic2001|FE1+FE2+YEARAPAY|0|School,temp[grepl("(Asian|Black|Hispanic|White)$",temp$EthGen),])))
    R<-tf(R,2+3*(type=="14"),summary(felm(PercApply~Asian1995+Asian1996+Asian1997+Asian1998+Asian1999+Asian2000+Asian2001+
                                            Black1995+Black1996+Black1997+Black1998+Black1999+Black2000+Black2001+
                                            Hispanic1995+Hispanic1996+Hispanic1997+Hispanic1998+Hispanic1999+Hispanic2000+Hispanic2001|FE1+FE2+YEARAPAY|0|School,temp[grepl("(Asian|Black|Hispanic|White)$",temp$EthGen),],weights = temp$Grads[grepl("(Asian|Black|Hispanic|White)$",temp$EthGen)])))
    R<-tf(R,3+3*(type=="14"),summary(felm(PercApply~Asian1995+Asian1996+Asian1997+Asian1998+Asian1999+Asian2000+Asian2001+
                                            Black1995+Black1996+Black1997+Black1998+Black1999+Black2000+Black2001+
                                            Hispanic1995+Hispanic1996+Hispanic1997+Hispanic1998+Hispanic1999+Hispanic2000+Hispanic2001|FE1+FE2+YEARAPAY|0|School,temp[grepl("(Asian|Black|Hispanic|White)_[MF]",temp$EthGen),],weights = temp$Grads[grepl("(Asian|Black|Hispanic|White)_[MF]",temp$EthGen)])))
  }
  write.csv(R,file="Figures/AA_ChangeApps_HS.csv")
}










###
#  Use IPUMS ACS data to see what's happening big-picture
###
if(IPUMS_Data){
  load(paste0(secure_derived,"AffAct_Data.Rda"))
    A<-A[A$CQUARTERAP==2&A$SCHTYPE%in%c("A","B")&!is.na(A$hsgpa),]
  load("Figures/Tables.Rda")
  
  load("Data/Derived/ACS_Data_Clean.Rda")
  
  #Start by estimating statistics: number of missing high-earning URM workers
  YoungAll<-data[data$AGE%in%30:34&data$YEAR==2014,]
  Young<-data[data$AGE%in%30:34&data$YEAR==2014&data$STATEFIP==6,]
  Young07<-data[data$AGE%in%30:34&data$YEAR==2007&data$STATEFIP==6,]
  
  aa<-sum(A$YEARAPAY%in%1998:2002&A$URM==1)*mean(as.numeric(gsub("[^0-9.-]","",R3[81,paste0("wage_sum_",12:16,"_100000")])))/100 #Number of additional rich URM youths, relative to 1995
  u<-sum(Young$PERWT[Young$URM==1&Young$INCWAGE>=100000])
  aa/(u-aa) ; aa/u
  weighted.quantile(YoungAll$INCWAGE[YoungAll$INCWAGE>0],YoungAll$PERWT[YoungAll$INCWAGE>0],.95)
  
  aa<-sum(A$YEARAPAY%in%1998:2002&A$URM==1)*mean(as.numeric(gsub("[^0-9.-]","",R3[81,paste0("wage_sum_",12:16,"_150000")])))/100 #Number of additional rich URM youths, relative to 1995
  u<-sum(Young$PERWT[Young$URM==1&Young$INCWAGE>=150000])
  aa/(u-aa) ; aa/u
  
  num<-sum(Young$PERWT[Young$URM==1])
  p<-u/num
  p07<-sum(Young07$PERWT[Young07$URM==1&Young07$INCWAGE>=100000])/sum(Young07$PERWT[Young07$URM==1])
  aa/num/p07 #(u-aa)/num/p07-p/p07 #This is the movement on the graph
  
  #Total Change in wages
  (sum(A$YEARAPAY%in%1998:2002&A$URM==1)*mean(as.numeric(gsub("[^0-9.-]","",R3[81,paste0("wage_sum_",12:16)]))))/
    sum(Young$Wage[Young$URM]*Young$PERWT[Young$URM]) #Loss of 0.5 p.p. of wages
  
  #Shares of all workers
  mean(data$URM[data$AGE%in%30:34&data$YEAR==2010&data$CA])
  mean(data$URM[data$AGE%in%30:34&data$YEAR==2010&data$CA&data$K100])
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces Figure A-13.
  ##########################################################################################################
  ##########################################################################################################
  for(shift in c(0,3)){
    for(i in c(25,30,35,40)) data$AgeGroup[data$AGE%in%((i+2*(i==35)+shift):(i+4-2*(i==25)+shift))]<-i
    for(c in c("CA","FromCA")){
      data$Cal<-data[,c]
      coldef<-64
      temp<-ddply(data[!is.na(data$AgeGroup)&data$YEAR%in%2001:2018&!is.na(data$Eth)&!is.na(data$Cal),],.(Eth,YEAR,Cal,AgeGroup),function(x){
        return(data.frame(K100=weighted.mean(x$K100,x$PERWT),K150=weighted.mean(x$K150,x$PERWT),GradSch=weighted.mean(x$GradSch,x$PERWT),K75=weighted.mean(x$INCWAGE>75000,x$PERWT),MedInc=weighted.median(x$INCWAGE,x$PERWT),MedInc_Condit=weighted.median(x$INCWAGE[x$INCWAGE>0],x$PERWT[x$INCWAGE>0]),STEM=weighted.mean(x$STEM),STEM_Condit=weighted.mean(x$STEM[x$EDUCD>=101],x$PERWT[x$EDUCD>=101]),MedInc_Coll=weighted.median(x$INCWAGE[x$EDUCD>coldef],x$PERWT[x$EDUCD>coldef]),K100_Coll=weighted.mean(x$K100[x$EDUCD>coldef],x$PERWT[x$EDUCD>coldef]),MedInc_NonColl=weighted.median(x$INCWAGE[inrange(x$EDUCD,65,100)],x$PERWT[inrange(x$EDUCD,65,100)]),K100_NonColl=weighted.mean(x$K100[inrange(x$EDUCD,65,100)],x$PERWT[inrange(x$EDUCD,65,100)])))
      }) #Gets noisy at 125k and above
      for(v in names(temp)[-(1:4)]) temp[,v]<-(temp[,v]+c(NA,temp[-nrow(temp),v]))/2 ; temp<-temp[temp$YEAR!=2001,] #Two-year moving average
      for(v in names(temp)[grepl("STEM",names(temp))]) temp[temp$YEAR==2007,v]<-temp[temp$YEAR==2009,v] #Works as long as the order stays basically the same... STEM doesn't appear in survey until 2009
      Baseline<-temp[temp$YEAR==2007+shift,!names(temp)=="YEAR"] ; names(Baseline)[-(1:3)]<-paste0(names(Baseline)[-(1:3)],"_B")
      for(v in names(Baseline)[-(1:3)]) Baseline[Baseline[,v]==0,v]<-NA
      t<-names(temp)[-(1:4)] ; temp<-merge(temp,Baseline) ; for(v in t) temp[,paste0(v,"_N")]<-temp[,v]/temp[,paste0(v,"_B")]
      temp<-temp[order(temp$Cal,temp$AgeGroup,temp$Eth,temp$YEAR),]
      for(v in names(temp)[grepl("STEM",names(temp))]) temp[temp$YEAR<2009,v]<-NA #Delete STEM info
      par(mar = c(5,4,2,2))
      
      if(c=="FromCA") fromCA<-temp
      if(c=="CA") CA<-temp
      vars<-list(list("MedInc_N",c(.6,1.2),paste0("% Rel. to ",2007+shift),.85),
                 list("MedInc_Coll_N",c(.6,1.2),paste0("% Rel. to ",2007+shift),.85),
                 list("MedInc_NonColl_N",c(.6,1.2),paste0("% Rel. to ",2007+shift),.85),
                 list("MedInc_Condit_N",c(.8,1.2),paste0("% Rel. to ",2007+shift),.85),
                 list("K100_N",c(.4,1.2)+.25*(shift==3),paste0("% Rel. to ",2007+shift),.72),
                 list("K100_Coll_N",c(.4,1.2)+.25*(shift==3),paste0("% Rel. to ",2007+shift),.72),
                 list("K100_NonColl_N",c(.4,1.2+.35*(shift==3))+.25*(shift==3),paste0("% Rel. to ",2007+shift),.72),
                 list("K75_N",c(.5,1.2),paste0("% Rel. to ",2007+shift),.85),
                 list("K150_N",c(.5,1.2),paste0("% Rel. to ",2007+shift),.85),
                 list("GradSch_N",c(.6,1.4),paste0("% Rel. to ",2007+shift),.85),
                 list("STEM_N",c(.9,1.6),paste0("% Rel. to ",2007+shift),.85),
                 list("STEM_Condit_N",c(.9,1.6),paste0("% Rel. to ",2007+shift),.85))
      for(v in vars){
        var<-v[[1]]
        extraname<-"" ; if(shift==3) extraname<-"_Shift3"
        tiff(paste0("Figures/ACS_",var[[1]],"_",c,extraname,".tiff"),5,3.5,"in",res=300)
          plot(2002:2018,temp[temp$Eth=="URM"&temp$Cal&temp$AgeGroup==30,var],type="l",ylim=v[[2]],xlab="ACS Response Year",ylab=v[[3]],lwd=3,cex.axis=1.5,cex.lab=1.5,mgp=c(2.5,1,0))
          polygon(c(2009.5,2009.5,2014.5,2014.5)+shift,c(0,100,100,0),col='gray80',lty = 0)
          polygon(c(2009.5,2009.5,2007.5,2007.5)+shift,c(0,100,100,0),col='gray90',lty = 0)
          abline(h=1,col="gray20")
          par(new = T)
          plot(2002:2018,temp[temp$Eth=="URM"&temp$Cal&temp$AgeGroup==30,var],type="l",ylim=v[[2]],xlab="ACS Response Year",ylab=v[[3]],lwd=3,cex.axis=1.5,cex.lab=1.5,mgp=c(2.5,1,0))
          lines(2002:2018,temp[temp$Eth=="URM"&temp$Cal&temp$AgeGroup==35,var],col=t_col("blue",50),lty=2,lwd=2)
          lines(2002:2018,temp[temp$Eth=="URM"&!temp$Cal&temp$AgeGroup==30,var],col=t_col("forestgreen",50),lty=1,lwd=2)
          lines(2002:2018,temp[temp$Eth=="White"&temp$Cal&temp$AgeGroup==30,var],col=t_col("firebrick2",50),lty=1,lwd=2)
          lines(2002:2018,temp[temp$Eth=="Asian"&temp$Cal&temp$AgeGroup==30,var],col=t_col("firebrick2",50),lty=2,lwd=2)
          text(x=2011+shift,y=v[[2]][2]+0.05,labels=paste0("Prop. 209 Phase-In\nfor ",30+shift,"-",34+shift," Cohorts"),cex = 1.15,pos=1)
        dev.off()
      }
    }
  }
  #Now just produce a legend figure
  tiff(paste0("Figures/ACS_Legend.tiff"),width=8,height=3,units="in",res=500)
    plot(100000,100000,bty="n",box=F,yaxt="n",xaxt="n",xlim=c(0,50),ylim=c(0,10),xlab="",ylab="")
    legend(1,8,legend=c("30-34 CA URM","37-39 CA URM","30-34 Non-CA URM","30-34 CA White","30-34 CA Asian"),lty=c(1,2,1,1,2),lwd=c(3,2,2,2,2),col=c("black",t_col("blue",50),t_col("forestgreen",50),t_col("firebrick2",50),t_col("firebrick2",50)),cex=.6,ncol = 5)
  dev.off()
  tiff(paste0("Figures/ACS_Legend_Shift3.tiff"),width=8,height=3,units="in",res=500)
    plot(100000,100000,bty="n",box=F,yaxt="n",xaxt="n",xlim=c(0,50),ylim=c(0,10),xlab="",ylab="")
    legend(1,8,legend=c("33-37 CA URM","40-42 CA URM","33-37 Non-CA URM","33-37 CA White","33-37 CA Asian"),lty=c(1,2,1,1,2),lwd=c(3,2,2,2,2),col=c("black",t_col("blue",50),t_col("forestgreen",50),t_col("firebrick2",50),t_col("firebrick2",50)),cex=.6,ncol = 5)
  dev.off()
}











if(CK_Replication){
  load(paste0(secure_cb,"SATUniverse.Rda"))
    load(paste0(secure_cb,"SAT_IDs.Rda"))
  load(paste0(secure_derived,"uadm_data_extract_processed.Rda"))
    data<-merge(Data,SAT_IDs) ; Data<-merge(Data,SAT_IDs,all.x=T)
  
  S<-S[S$Year%in%1994:2001,]
  S$URM<-S$Ethnicity%in%c("Black","Hispanic","Nat_Amer")
  S$URM_Pred<-S$Ethnicity_Pred%in%c("Black","Hispanic","Nat_Amer")
  S$SAT<-S$SATIV_CB+S$SATIM_CB
  S$SAT_Cat<-0 ; S$SAT_Cat[S$SAT>1150]<-1 ; S$SAT_Cat[S$SAT>1300]<-2
  S$CLASRANK<-S$Rank_cb ; S$CLASRANK[is.na(S$CLASRANK)]<-999
  S$hsgpa_censor<-S$hsgpa_cb ; S$hsgpa_censor[NA_to_F(S$hsgpa>4)]<-4
  S$AIS<-round((S$hsgpa_censor*1000+S$SATIV+S$SATIM+S$SATIIW_CB+S$SATIIM_CB+S$SATIIT_CB)/10)*10
  data$APPUC<-(data$APP01==1|data$APP03==1|data$APP04==1|data$APP05==1|data$APP06==1|data$APP07==1|data$APP08==1|data$APP09==1)
  S<-merge(S,data[,grepl("^APP[0-9U]|IDCB",names(data))],all.x=T)
  for(v in names(S)[grepl("^APP[0-9U]",names(S))]) S[is.na(S[,v]),v]<-0
  S$APP14<-(S$APP01==1|S$APP04==1) ; S$Send14<-(S$Send01==1|S$Send04==1)
  S$URM678<-S$URM*(S$Year%in%c(1996,1997,1998)) ; S$URM901<-S$URM*(S$Year%in%c(1999,2000,2001))
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces the statistics presented in Tables F-1, F-2, F-3, and F-4.
  ##########################################################################################################
  ##########################################################################################################
  R<-data.frame(matrix("",100,25),stringsAsFactors = F)
  tf<-function(R,c,reg,reg1=NULL){
    if(sum(grepl("URM_Pred",row.names(reg$coefficients)))>0) urm<-"URM_Pred" else urm<-"URM"
    r<-reg$coefficients[paste0(urm,"TRUE:factor(Year)",1995:2001),]
    R[seq(1,19,3),c]<-ex(r[,1],3)
    R[seq(2,20,3),c]<-ex(r[,2],3,par=T)
    
    R[22,c]<-"X"
    if(urm=="URM_Pred") R[23,c]<-"X"
    
    if(!is.null(reg1)){
      R[25,c]<-ex(reg1$coefficients[row.names(reg1$coefficients)=="URM901",1],3)
      R[26,c]<-ex(reg1$coefficients[row.names(reg1$coefficients)=="URM901",2],3,par=T)
    }
    
    
    R[28,c]<-ex(reg$r2,2,comma=T)
    R[29,c]<-ex(reg$N,0,comma=T)
    return(R)
  }
  R<-tf(R,1,summary(felm(SendUC~URM*factor(Year)|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S,exactDOF = T),robust=T))
  R<-tf(R,2,summary(felm(APPUC~URM*factor(Year)|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S,exactDOF = T),robust=T))
  R<-tf(R,3,summary(felm(SendUC~URM*factor(Year)|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$SAT>1150,],exactDOF = T),robust=T))
  R<-tf(R,4,summary(felm(APPUC~URM*factor(Year)|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$SAT>1150,],exactDOF = T),robust=T))
  R<-tf(R,5,summary(felm(SendUC~URM*factor(Year)|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$hsgpa_cb>3.99,],exactDOF = T),robust=T),summary(felm(SendUC~factor(Year)+URM678+URM901|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$hsgpa_cb>3.99,],exactDOF = T),robust=T))
  R<-tf(R,6,summary(felm(APPUC~URM*factor(Year)|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$hsgpa_cb>3.99,],exactDOF = T),robust=T),summary(felm(APPUC~factor(Year)+URM678+URM901|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$hsgpa_cb>3.99,],exactDOF = T),robust=T))
  R<-tf(R,7,summary(felm(APPUC~URM_Pred*factor(Year)|Ethnicity_Pred+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$hsgpa_cb>3.99,],exactDOF = T),robust=T))
  
  R<-tf(R,9,summary(felm(Send14~URM*factor(Year)|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S,exactDOF = T),robust=T))
  R<-tf(R,10,summary(felm(APP14~URM*factor(Year)|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S,exactDOF = T),robust=T))
  R<-tf(R,11,summary(felm(Send14~URM*factor(Year)|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$SAT>1150,],exactDOF = T),robust=T))
  R<-tf(R,12,summary(felm(APP14~URM*factor(Year)|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$SAT>1150,],exactDOF = T),robust=T))
  R<-tf(R,13,summary(felm(Send14~URM*factor(Year)|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$hsgpa_cb>3.99,],exactDOF = T),robust=T),summary(felm(Send14~factor(Year)+URM678+URM901|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$hsgpa_cb>3.99,],exactDOF = T),robust=T))
  R<-tf(R,14,summary(felm(APP14~URM*factor(Year)|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$hsgpa_cb>3.99,],exactDOF = T),robust=T),summary(felm(APP14~factor(Year)+URM678+URM901|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$hsgpa_cb>3.99,],exactDOF = T),robust=T))
  R<-tf(R,15,summary(felm(Send14~URM*factor(Year)|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$AIS%in%5500:7000,],exactDOF = T),robust=T))
  R<-tf(R,16,summary(felm(APP14~URM*factor(Year)|Ethnicity+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$AIS%in%5500:7000,],exactDOF = T),robust=T))
  R<-tf(R,17,summary(felm(APP14~URM_Pred*factor(Year)|Ethnicity_Pred+Year+SAT_Cat+FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$AIS%in%5500:7000,],exactDOF = T),robust=T))
  S$SAT_Norm<-(S$SAT-mean(S$SAT,na.rm=T))/sd(S$SAT,na.rm=T)
  hi<-summary(felm(APP14~Send14*URM*factor(Year)*SAT_Norm|FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S,exactDOF = T),robust=T)
  R[1:63,18]<-row.names(hi$coefficients)
  R[1:63,19]<-ex(hi$coefficients[,1],3)
  R[1:63,20]<-ex(hi$coefficients[,2],3,par=T)
  R[1:63,21]<-ex(hi$coefficients[,4],3)
  R[65,19]<-ex(hi$r.squared,2)
  R[66,19]<-ex(hi$N,0,comma=T)
  write.csv(R,file="Figures/CK_Regs.csv")
  
  #Summary statistics
  cor(S$Send01,S$APP01)
  table(S$Send01,S$APP01)
  summary(felm(APP14~Send14*URM*factor(Year)*SAT_Norm|FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S,exactDOF = T),robust=T)
  summary(felm(APP14~URM*factor(Year)*SAT_Norm|FATHEDUC+MOTHEDUC+hsgpa_cb+CLASRANK+Grade_cb,data=S[S$Send14,],exactDOF = T),robust=T)
}









if(Value_Added_Coefficients){
  load(paste0(secure_nsc,"Six_Year_Grad/Full_NSC.Rda"))
  load(paste0(secure_derived,"NSC_ValueAdded_Measures.Rda"))
    names(FEs)[1]<-"OPEID"
    
  a<-FEs[!is.na(FEs$OPEID),!grepl("SEX|CETH|Eth|College_Name|AIS_Bin",names(FEs))]
    FEs<-FEs[!is.na(FEs$NSC_EnrFirstIncCC_College_Code_Eth),]
    FEs$OPEID<-gsub("\\s-\\s.*$","",FEs$NSC_EnrFirstIncCC_College_Code_Eth)
    fe<-FEs[FEs$CETHNICA_Cat=="Black",grepl("OPEID|FE_[WU]",names(FEs))]
    names(fe)[-1]<-paste0(names(fe)[-1],"_Black")
  a<-merge(a,fe,all.x=T)
    fe<-FEs[FEs$CETHNICA_Cat=="Hispanic",grepl("OPEID|FE_[WU]",names(FEs))]
    names(fe)[-1]<-paste0(names(fe)[-1],"_Hispanic")
  a<-merge(a,fe,all.x=T)  
  
  nsc<-read.csv("Data/Raw/IPEDS_NSC_Compare.csv",stringsAsFactors = F)
    nsc$OPEID<-gsub("^0*|..\\s*$","",nsc$Office.of.Postsecondary.Education..OPE..ID.Number..HD2008.)
    nsc<-nsc[NA_to_F(nsc$OPEID!=""),c("Institution.Name","OPEID")]
    nsc<-nsc[!(allduplicated(nsc$OPEID)&grepl("District|Office",nsc$Institution.Name)),] #Wrong names
    nsc<-ddply(nsc,.(OPEID),function(x) data.frame(Institution.Name=x$Institution.Name[1]))
  a<-merge(a,nsc,all=T)
  names(a)<-c("OPEID",names(a)[grepl("FE",names(a))],"InstName")
  a<-a[order(a$FE_GradVA_Chetty_AA,decreasing=T),c("InstName",names(a)[grepl("FE",names(a))],"OPEID")]
  
  #Re-norm FEs to be relative to CSU Long Beach
  for(v in names(a)[grepl("FE.*(Mountjoy|Chetty|DKCB|Raw)(_AA)?(_Black|_Hispanic)?$",names(a))&!grepl("se|obs",names(a))]) a[,v]<-a[,v]-a[NA_to_F(a$InstName=="California State University-Long Beach"),v]
  for(v in names(a)[grepl("(FE.*GradVA|^FE_GradRate)",names(a))&!grepl("obs",names(a))]) a[,v]<-a[,v]*100
  
  for(v in c("FE_GradVA_Mountjoy_AA","FE_GradVA_Raw_AA","FE_GradVA_DKCB_AA","FE_GradVA_Chetty_AA")){ ; a[,v]<-ex(a[,v],1) ; a[grepl("N",a[,v]),v]<-'="-"' ; }
  for(v in c("FE_WageVA_Mountjoy_AA","FE_WageVA_Raw_AA","FE_WageVA_DKCB_AA","FE_WageVA_Chetty_AA","FE_WageVA_Chetty_AA_Black","FE_WageVA_Chetty_AA_Hispanic")){ ; a[,v]<-ex(round(a[,v]/100)*100,0,comma=T) ; a[grepl("N",a[,v]),v]<-'="-"' ; }
  for(v in c("FE_hsgpa_Mountjoy_AA","FE_hsgpa_Raw_AA","FE_hsgpa_Chetty_AA")){ ; a[,v]<-ex(a[,v],2,comma=T) ; a[grepl("N",a[,v]),v]<-'="-"' ; }  
  for(v in c("FE_obs_GradVA_Chetty_AA","FE_obs_WageVA_Chetty_AA")){ ; a[,v]<-ex(a[,v],0,comma=T) ; a[grepl("N",a[,v]),v]<-'="-"' ; }
  a$InstName<-as.character(a$InstName)
  a<-a[order(a$InstName),]
  
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces the statistics presented in Tables I-1, I-2, and I-3.
  ##########################################################################################################
  ##########################################################################################################
  vars<-list(c("InstName","FE_GradVA_Raw_AA","FE_GradVA_Mountjoy_AA","FE_GradVA_Chetty_AA","FE_WageVA_Raw_AA","FE_WageVA_Mountjoy_AA","FE_WageVA_Chetty_AA","FE_WageVA_Chetty_AA_Black","FE_WageVA_Chetty_AA_Hispanic","FE_hsgpa_Raw_AA","FE_hsgpa_Mountjoy_AA","FE_hsgpa_Chetty_AA","FE_obs_WageVA_Chetty_AA"))
  for(v in 1:length(vars)){
    uc<-a[grepl("^University of Calif",a$InstName)&!is.na(a$FE_GradVA_Raw_AA),vars[[v]]]
    uc$InstName<-gsub("University of California-","UC ",uc$InstName) ; uc<-uc[uc$FE_GradVA_Chetty_AA!='=""',] ; uc<-uc[order(uc$InstName),] ; uc<-cbind(uc[1:4,],uc[5:8,])
    
    csu<-a[grepl("CALIFORNIA STATE UNIV|CALIFORNIA POLYTECHNIC STATE|CALIFORNIA STATE POLY|SAN DIEGO STATE|SAN FRANCISCO STATE UN|SAN JOSE STATE UNIV|HUMBOLDT STATE UNI|SONOMA STATE UNIV|CALIFORNIA MARITIME ACADEMY",a$InstName,ignore.case=T),vars[[v]]]
    csu<-csu[csu$FE_GradVA_Chetty_AA!='=""',] ; csu<-csu[order(csu$InstName),] ; csu<-cbind(csu[1:10,],csu[11:20,])
    
    cc1<-unique(Full_NSC$`College.Code/Branch`[Full_NSC$`Public./.Private`=="Public"&Full_NSC$`College.State`=="CA"&Full_NSC$`2-year./.4-year`==2])
    cc1<-gsub("^0*|...\\s*$","",cc1) ; cc1<-cc1[!duplicated(cc1)]
    cc<-a[a$OPEID%in%cc1,vars[[v]]]
    cc<-cc[grepl("[0-9]",cc$FE_GradVA_Chetty_AA),]
    cc<-cc[order(cc$InstName),] ; cc<-cbind(cc[1:25,],cc[26:50,])
    
    priv<-a[!(grepl("^University of Calif",a$InstName)&!is.na(a$FE_GradVA_Raw_AA))&
              !grepl("CALIFORNIA STATE UNIV|CALIFORNIA POLYTECHNIC STATE|CALIFORNIA STATE POLY|SAN DIEGO STATE|SAN FRANCISCO STATE UN|SAN JOSE STATE UNIV|HUMBOLDT STATE UNI|SONOMA STATE UNIV|CALIFORNIA MARITIME ACADEMY",a$InstName,ignore.case=T)&
              !(a$OPEID%in%cc1)&
              grepl("[0-9]",a$FE_GradVA_Raw_AA),vars[[v]]]
    priv<-priv[!priv$FE_GradVA_Chetty_AA=='=""',] ; priv<-priv[order(priv$InstName),] ; priv<-cbind(priv[1:(nrow(priv)/2),],priv[(nrow(priv)/2+1):nrow(priv),])
    
    R<-rbind(uc,csu,cc,priv)
    for(i in 1:ncol(R)) R[is.na(R[,i]),i]<-'=""'
    for(i in c(1,6)){
      R[,i]<-gsub("University of California-","UC ",R[,i])
      R[,i]<-gsub("California State University-","CSU ",R[,i])
      R[,i]<-gsub("California","CA",R[,i])
      R[,i]<-gsub("Polytechnic","Poly.",R[,i])
      R[,i]<-gsub("Technical","Tech.",R[,i])
      R[,i]<-gsub("University","Univ.",R[,i])
      R[,i]<-gsub("Community College","C.C.",R[,i])
      R[,i]<-gsub("College","C.",R[,i])
      R[,i]<-gsub("Los Angeles","LA",R[,i])
      R[,i]<-gsub("\\s*District|[- ]*Administrative Office|(-[^0-9]|Office).*$","",R[,i])
    }
    write.csv(R,file="Figures/NSC_GradRateList.csv")
  }
}






if(Miscellaneous){
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces the statistics presented in Tables A-1 and A-2.
  ##########################################################################################################
  ##########################################################################################################
  load(paste0(secure_derived,"AffAct_Data.Rda"))
    A<-A[A$CQUARTERAP==2&A$SCHTYPE%in%c("A","B")&!is.na(A$hsgpa)&A$YEARAPAY%in%1994:2002,]
  hi<-data.frame(c(A$NSC_Major1,A$NSC_Major2,A$NSC_Major3),1) ; names(hi)<-c("Majors","One")
    hi<-hi[!hi[,1]%in%c("","NA","NOT APPLICABLE"),]
    hi<-aggregate(One~Majors,hi,sum)
  load("Data/Derived/NSC_STEM_fromCIP.Rda")
    hi$STEM<-hi$Majors%in%STEM
    hi<-hi[order(hi$STEM,hi$One,decreasing = T),]
    hi<-cbind(hi[hi$STEM,c("Majors","One")][1:50,],hi[hi$STEM,c("Majors","One")][51:100,],hi[!hi$STEM,c("Majors","One")][1:50,],hi[!hi$STEM,c("Majors","One")][51:100,])
  write.csv(hi,file="Figures/STEM_List.csv")
  
  
  
  
  
  
  #Take a look at aggregate university enrollment
  ipeds<-data.frame(read_csv("Data/Raw/IPEDS_Agg_Enroll.csv"))
  ipeds$OPEID<-as.integer(gsub("..$","",ipeds$Office.of.Postsecondary.Education.ID.Number..FA2000HD.))
  ipeds$Level<-ipeds$Level.of.institution..HD2002.
  ipeds$Sector<-ipeds$Control.of.institution..HD2002.
  ipeds$State<-ipeds$USPS.state.abbreviation..HD2002.
  ipeds$GradRate<-ipeds$Grand.total..GR1997..4.year.institutions..Completers.within.150..of.normal.time./ipeds$Grand.total..GR1997..4.year.institutions..Adjusted.cohort..revised.cohort.minus.exclusions..*100
  temp<-data.frame()
  for(i in 1994:2002){ #Reshape long
    hi<-ipeds[,grepl(paste0("UnitID|Institution|OPEID|Sector|Level|State|GradRate|",i,"|EF",i-1900),names(ipeds))]
    hi$Year<-i
    hi$Total<-hi[,grepl("Grand",names(hi))&!grepl("institutions",names(hi))]
    hi$URM<-hi[,grepl("Black",names(hi))]+hi[,grepl("^Hispanic",names(hi))]+hi[,grepl("Indian",names(hi))]
    hi$NonURM<-hi$Total-hi$URM
    hi<-hi[,!grepl("[0-9]",names(hi))]
    temp<-rbind(temp,hi)
  } ; ipeds<-temp
  ipeds<-ipeds[NA_to_F(ipeds$Total>0),]
  
  #Merge in value-added estimates
  load(paste0(secure_derived,"NSC_ValueAdded_Measures.Rda"))
    FEs$OPEID<-FEs$NSC_EnrFirstIncCC_College_Code
    FEs<-FEs[FEs$OPEID%in%ipeds$OPEID&!is.na(FEs$OPEID),]
    ipeds<-merge(ipeds,FEs[,c("OPEID","FE_WageVA_Chetty_AA")],all.x=T)
    
  ##########################################################################################################
  ##########################################################################################################
  #The following code produces Figure A-2.
  ##########################################################################################################
  ##########################################################################################################
  hi<-aggregate(.~Year,ipeds[ipeds$State=="CA"&ipeds$Sector==1&ipeds$Level==1,c("Year","Total","URM","NonURM")],sum) 
  for(v in c("Total","URM","NonURM")) hi[,v]<-(hi[,v]-c(NA,hi[-nrow(hi),v]))/hi[,v]*100
  tiff("Figures/UnivEnroll_CAPubFourYears_NumStud.tiff",5,3.5,units="in",res=300)
    plot(hi$Year,hi$Total,ylim=c(-10,12),pch=16,cex.lab=1.5,cex.axis=1.2, xlab="", ylab="Δ Enrollment (%)",mgp=c(2.5,1,0))
    polygon(c(1997.5,1997.5,1998.5,1998.5),c(-100,100,100,-100),col='gray90',lty = 0)
    abline(h=0,col="gray50",lty=3)
    box()
    points(hi$Year,hi$Total,pch=16)
    points(hi$Year,hi$URM,pch=7)
    points(hi$Year,hi$NonURM,lty=3,pch=2)
    legend(2000.2,-1.1,legend=c("All","URM","Non-URM"),pch=c(16,7,2),cex=.8,ncol = 1)
  dev.off()
  
  #CA all four-year enrollment
  hi<-aggregate(.~Year,ipeds[ipeds$State=="CA"&ipeds$Level==1,c("Year","Total","URM","NonURM")],sum) 
  for(v in c("Total","URM","NonURM")) hi[,v]<-(hi[,v]-c(NA,hi[-nrow(hi),v]))/hi[,v]*100
  tiff("Figures/UnivEnroll_CAFourYears_NumStud.tiff",5,3.5,units="in",res=300)
    plot(hi$Year,hi$Total,ylim=c(-10,12),pch=16,cex.lab=1.5,cex.axis=1.2, xlab="", ylab="Δ Enrollment (%)",mgp=c(2.5,1,0))
    polygon(c(1997.5,1997.5,1998.5,1998.5),c(-100,100,100,-100),col='gray90',lty = 0)
    abline(h=0,col="gray50",lty=3)
    box()
    points(hi$Year,hi$Total,pch=16)
    points(hi$Year,hi$URM,pch=7)
    points(hi$Year,hi$NonURM,lty=3,pch=2)
    legend(2000.2,-1.1,legend=c("All","URM","Non-URM"),pch=c(16,7,2),cex=.8,ncol = 1)
  dev.off()
    
  #Weighted avg VA, CA four-year publics
  for(v in c("Total","URM","NonURM")) ipeds[,paste0(v,"_VA")]<-ipeds[,v]*ipeds$FE_WageVA_Chetty_AA
  hi<-aggregate(.~Year,ipeds[ipeds$State=="CA"&ipeds$Sector==1&ipeds$Level==1,c("Year","Total_VA","URM_VA","NonURM_VA","Total","URM","NonURM")],sum)
  for(v in c("Total_VA","URM_VA","NonURM_VA")){
    hi[,v]<-hi[,v]/hi[,gsub("...$","",v)]
    hi[,v]<-(hi[,v]-c(NA,hi[-nrow(hi),v])) #/hi[,v]*100
  }
  tiff("Figures/UnivEnroll_CAPubFourYears_AvgVA.tiff",5,3.5,units="in",res=300)
    plot(hi$Year,hi$Total_VA/1000,ylim=c(-1,.5),pch=16,cex.lab=1.5,cex.axis=1.2, xlab="", ylab="Δ Dollars ('000s)",mgp=c(2.5,1,0))
    polygon(c(1997.5,1997.5,1998.5,1998.5),c(-100,100,100,-100),col='gray90',lty = 0)
    abline(h=0,col="gray50",lty=3)
    box()
    points(hi$Year,hi$Total_VA/1000,pch=16)
    points(hi$Year,hi$URM_VA/1000,pch=7)
    points(hi$Year,hi$NonURM_VA/1000,lty=3,pch=2)
    legend(2000.2,-.400,legend=c("All","URM","Non-URM"),pch=c(16,7,2),cex=.8,ncol = 1)
  dev.off()
  
  #Weighted avg VA, CA four-year all
  for(v in c("Total","URM","NonURM")) ipeds[,paste0(v,"_VA")]<-ipeds[,v]*ipeds$FE_WageVA_Chetty_AA
  hi<-aggregate(.~Year,ipeds[ipeds$State=="CA"&ipeds$Level==1,c("Year","Total_VA","URM_VA","NonURM_VA","Total","URM","NonURM")],sum)
  for(v in c("Total_VA","URM_VA","NonURM_VA")){
    hi[,v]<-hi[,v]/hi[,gsub("...$","",v)]
    hi[,v]<-(hi[,v]-c(NA,hi[-nrow(hi),v])) #/hi[,v]*100
  } 
  tiff("Figures/UnivEnroll_CAFourYears_AvgVA.tiff",5,3.5,units="in",res=300)
    plot(hi$Year,hi$Total_VA/1000,ylim=c(-1,.5),pch=16,cex.lab=1.5,cex.axis=1.2, xlab="", ylab="Δ Dollars ('000s)",mgp=c(2.5,1,0))
    polygon(c(1997.5,1997.5,1998.5,1998.5),c(-100,100,100,-100),col='gray90',lty = 0)
    abline(h=0,col="gray50",lty=3)
    box()
    points(hi$Year,hi$Total_VA/1000,pch=16)
    points(hi$Year,hi$URM_VA/1000,pch=7)
    points(hi$Year,hi$NonURM_VA/1000,lty=3,pch=2)
    legend(2000.2,-.4,legend=c("All","URM","Non-URM"),pch=c(16,7,2),cex=.8,ncol = 1)
  dev.off()
  
  ##########################################################################################################
  ##########################################################################################################
  #Figure A-3 is archival, and is not produced here.
  ##########################################################################################################
  ##########################################################################################################
}
